2.2.4. Thermal Range of High Photosynthesis (ΔT90)

The ΔT90 is the extent of the temperature range where assimilation rates are ≥90% of their seasonal maximum value calculated above, Anet,max. Considering the previously computed Anet,max and, at each time step, the average Anet over the whole canopy, we identified the maximum and minimum TL for which Anet is equal to 90% of Anet,max. ΔT90 corresponds to the difference between those maximum and minimum TL. We interpreted ΔT90 as a proxy of leaf thermotolerance: broad ΔT90 ranges mean that the CO2 assimilation rate is close to its maximum under a wide range of thermal conditions, while narrow ΔT90 range necessarily implies a higher likelihood to experience sub-optimal thermal conditions for assimilation.

#### *2.3. Plant Traits Selection and Determination of Their Role*

#### 2.3.1. Selection of Plant Traits

Only traits expected to have a direct most notable effect on leaf thermoregulation and C-fixation in APES were included in the analyses. Under set environmental conditions, TL is strongly affected by leaf size and shape [8], albedo [9] and traits related to stomatal conductance and hence regulation of transpiration [10]. These traits also affect the net CO2 assimilation rate, directly via its maximum rate, and indirectly, via the stomatal opening and the effects of temperature on leaf metabolic rates. Consequently, in the analyses below, we focused on the primary plant traits regulating the CO2 assimilation rate, the boundary layer and stomatal conductances (gb and gs, respectively), the absorption of solar radiation.

A total of six traits were selected: the maximum carboxylation rate at 25 ◦C, VCMAX,25; the two parameters of the stomatal model, i.e., stomatal model slope in well-watered conditions g1 and a parameter describing the sensitivity of g1 to soil water potential β [33]; the effective leaf thickness, lt, and PAR and NIR albedos, αPAR and αNIR. This selection was based on the following considerations. Assuming leaves as flat plates, the effective leaf thickness, lt, controls the boundary layer conductance gb [36,37]. Conversely, the stomatal conductance, gs, is mostly dependent on hydraulic and photosynthetic traits: in the stomatal model employed in APES [33], these aspects are summarized by the parameters g1 and β. Leaf temperature and photosynthesis also depend on the absorbed radiation, and hence on the leaf albedo for both photosynthetically active and near infrared radiation (αPAR and αNIR respectively). Finally, according to the Farquhar model, the net CO2 assimilation rate depends on the maximum carboxylation rate, VCMAX, the maximum electron transport rate, JMAX, and the day respiration rates, Rd. These rates are functions of the corresponding rates at 25 ◦C (VCMAX,25, JMAX,25 and Rd,25) and their response to temperature [31]. JMAX,25 and Rd,25 are often well correlated with VCMAX,25, following Medlyn et al. [38] and Launiainen et al. [29], respectively, thus allowing to consider only VCMAX,25 as a key parameter of photosynthesis. Further details on trait selection and their role in leaf temperature regulation are given in the Supplementary Materials.

For each of these six traits, we considered realistic ranges, based on a literature search. These ranges were maintained intentionally broad, beyond those typical of current dominant coniferous and deciduous tree species in boreal forests (Table 1). Conversely, the understory and forest floor features were fixed in the model. Concretely, the understory consists of seedlings of Norway spruce, Silver birch and other deciduous species (LAI <sup>∼</sup>0.5 m2 <sup>m</sup><sup>−</sup>2), as specified in [29].


**Table 1.** Summary of the traits considered in the General Sensitivity Analysis (GSA), their ranges of variation emerging from the literature globally and for mid-to-high latitudes and corresponding supporting references.

\* No specific values for mid-to-high latitude have been reported. However, boreal forests are dominated by species with needles. \*\* No specific values for mid-to-high latitude have been reported.

#### 2.3.2. General Sensitivity Analysis (GSA) and Determination of Dominant Traits

A General Sensitivity Analysis (GSA) was run on APES to identify which of the focal traits are most important to maximize both instantaneous and cumulative Anet and regulate TL, reducing the risk of high temperature damage. To this aim, 100 trait sets were randomly generated from uniform distributions of each trait over the ranges specified in Table 1. The traits were sampled independently. The model was run using each trait set for current growing conditions (Section 2.4.1). The obtained simulated time series of Anet and TL allowed the computation of the four metrics described in Section 2.2 for each trait set (maximum leaf temperature, TL,max; the thermal range of high photosynthesis, ΔT90, cumulated and maximum net assimilation, Anet,cum and Anet,max).

The model results were then grouped based on whether the metrics exceeded pre-set thresholds. Specifically, a threshold was set for each metric in order to split the 100 trait sets into two groups: one containing the m trait sets for which metric is above the threshold (Group 1) and the others *n* = (100 – *m*) (Group 2). Then, for each trait, the cumulative probability distribution for the *m* trait values in Group 1 and *n* trait values in Group 2 were obtained. The higher the distance between the cumulative distribution functions, the more important the role of the trait value in defining the outcome. To measure the distance between the two empirical functions, the Kolmogorov–Smirnov (KS) two-sample test was used. This nonparametric method allows comparing two samples and quantifying the distance between the empirical distributions of the two samples [51]; thus, it can be used to rank the traits according to how dominant they are in determining the final outcome [52–54].

The thresholds for the grouping were selected to be the highest ones that allowed the greatest number of trait sets in the Group 1 while the number of dominant traits remained unchanged [53]. Indeed, choosing an excessively low threshold leads to (almost) no distinction between Group 1 and Group 2 so that the separation between the two distribution functions is negligible for all traits (i.e., no trait emerges as dominant). Conversely, choosing an excessively high threshold reduces substantially the number of trait sets for which the metric is above the threshold: in this case, Group 1 might not be representative except for a minority of trait sets, leading to a substantial separation between the two distribution functions for all traits (i.e., all traits are considered dominant). This approach was followed for each metric except TL,max. For the latter, we used TCRIT—i.e., the temperature at which first symptoms of thermal damage appears. Following O'Sullivan et al. [3], based on the latitude of the reference study site (Section 2.4.1), TCRIT was set at 44.0 ◦C. This threshold also allowed a balanced partitioning of the 100 parameter sets into Group 1 and Group 2, i.e., the approach delineated above

would have led to a similar value for the threshold. The resulting thresholds were: TL,max = 44.0 ◦C, Anet,cum = 1,274 μmol m<sup>−</sup>2, Anet,max = 0.466 μmol m<sup>−</sup>2s−<sup>1</sup> and ΔT90 =6.3 ◦C.

#### 2.3.3. Selection of Focal Trait Combinations

Once the dominant traits were identified via Kolmogorov–Smirnov test, four sets of parameters were chosen to represent widely different outcomes in terms of thermotolerance and C fixation (Table 2). Given the focus on boreal forests, the focal trait combinations were chosen among those including trait values within the ranges typically observed in mid-to-high latitude forests (Table 1). We selected four trait combinations that showed contrasting thermal and productivity responses under identical growing conditions (2005 growing season). Specifically, trait combination 1 corresponds to the set of traits with the lowest performance (lowest Anet,cum and TL,MAX exceeding TCRIT in 2005): it is representative of those trait combinations with low productivity and high risk of thermal damage, at least under the current conditions (Low Productivity-High Risk; hereafter, LP-HR). Combination 2 was instead the best performing combination, leading to high assimilation rates while preventing thermal damage, i.e., it is representative of individuals with high productivity and low risk of thermal damage (High Productivity-Low Risk; HP-LR). Combination 3 led to high assimilation rates but also to TL,MAX on par with that of combination 1, thus representing trait sets leading to high productivity but also high risk of thermal damage (High Productivity-High Risk; HP-HR). Finally, combination 4 resulted in assimilation rates similar to combination 1 but it did prevent the occurrence of thermal damage; thus, it corresponds to low productivity but also low risk of thermal damage (Low Productivity-Low Risk). These differences in assimilation and thermal metrics are the result of the following differences between the trait values. Combinations 1 (LP-HR) and 4 (LP-LR) have similar and low VCMAX,25; combinations 2 (HP-LR) and 3 (HP-HR) have a similar VCMAX,25, but these are higher than combinations 1 and 4. Conversely, lt values in combinations 1 (LP-LR) and 3 (HP-HR) are similar and higher than the ones in combinations 2 (HP-LR) and 4 (LP-LR), the latter of which also have similar values. No patterns emerged relatively to the other traits, which were indeed not dominant as shown by the sensitivity analysis (see results in Section 3.1). These four combinations of traits together with their corresponding assimilation and thermal metrics are listed in Table 2.

**Table 2.** Combination of traits used to perform the climate scenarios analyses and the resulting values of the assimilation and thermal metrics for the growing season 2005. LP-HR: Low Productivity-High Risk of thermal damage; HP-LR: High Productivity-Low Risk of thermal damage; HP-HR: High Productivity-High Risk of thermal damage; LP-LR: Low Productivity-Low Risk of thermal damage.


## *2.4. Growing Conditions*

The GSA was performed under current climatic conditions, corresponding to selected years in a study site in Finland (Section 2.4.1). Conversely, the focal combinations of traits were tested under both current and future warmer and/or drier growing conditions, as defined in Section 2.4.2. This approach allowed us to determine whether thermal damage and suboptimal photosynthetic rates would occur given certain temperature increase and intensity of the water stress, thus allowing the quantification of potential effects of climate change and the importance of trait selection in the future.

#### 2.4.1. Current Growing Conditions

Current growing conditions are based on meteorological data observed in the Hyytiälä field station SMEAR II (Station for Measuring Forest Ecosystem-Atmosphere Relations), located in Hyytiälä, Southern Finland (61◦51' N, 24◦17' E). The annual mean temperature is 3 ◦C and total mean annual precipitation is 700 mm. All the meteorological data necessary for APES are routinely measured at this site. A complete description of the site and its measurements can be found in Kolari et al. [55]; Kulmala et al. [56] and Launiainen et al. [29].

The model was applied over selected growing seasons (2005 and 2006), here defined as the period of May 1st to September 30th [29], with the warmer period falling in the first half of July. The growing season 2005 can be considered as a normal one in the region, while 2006 was a dry year [57], slightly warmer than 2005. The two growing seasons were also characterized by a different precipitation timing: in 2005, the rainiest period occurred in late July and early August and the highest precipitation was recorded at the end of the growing season; in 2006, a dry spell occurred in late July and early August while two of the rainiest peaks were recorded in early July. Therefore, 2005 and 2006 present contrasting plant water availability, allowing us to analyze if these difference might lead to significant changes in thermoregulation and net CO2 assimilation.

#### 2.4.2. Scenarios of Future Growing Conditions

Boreal regions are experiencing a rapid warming and an increase in the frequency and intensity of periods with low water availability as a consequence of changes in the timing and magnitude of precipitation, as well as increases in evapotranspiration from longer and warmer growing seasons [15]. Therefore, to generate climate change scenarios, we focused on air temperature (TA) and soil water potential (ΨS).

The future scenarios were developed as follows. We considered two single days within the growing season of year 2005 as baseline: (i) the day where TL reached its maximum value (July 6th) and (ii) an average day with TA close to the median value for the growing season (May 20th). Then, we built scenarios of future growing conditions by shifting upwards the observed TA series of 1 to 10 ◦C combined with decreasing soil water potential ranging from −0.5 to −3 MPa. The soil water potential is considered constant within the day. For each scenario, we computed the instantaneous Anet and TL at midday. While acknowledging their importance, other environmental variables such as wind speed and relative humidity were not modified. Thus, the vapor pressure deficit (VPD) was computed considering the new temperature values and the unchanged relative humidity.

Furthermore, we analyzed how trait combinations might affect soil water availability and how such changes may alter maximum and cumulated net CO2 assimilation. To this aim, the shift in TA series by 1 to 10 ◦C was combined with imposing 5 to 30 day dry spells, with no precipitation. The dry spell was assumed to begin after the last precipitation event before the warmest period registered within the growing season of 2005 (i.e., June 27th). We then computed the corresponding Ψ<sup>s</sup> at the end of the 5–30 day dry period, and the Anet,cum over the dry down.

In all cases, the length of the growing season was maintained unaltered, even though global warming may lead to longer growing seasons in boreal regions [58,59]. This simplifying assumption bears no consequence of the results discussed below, because the aim was to contrast different trait combinations under the same abiotic conditions.

#### **3. Results**

#### *3.1. Dominant Traits Emerging from the General Sensitivity Analysis (GSA)*

The General Sensitivity Analysis (GSA) allowed the identification of the traits most relevant for the four metrics: TL,max, Anet,cum, Anet,max and ΔT90. The results are presented in Figure 1 and summarized next.

**Figure 1.** Results of the Global Sensitivity Analysis. Rows refer to the four metrics; columns to the six traits investigated. In each main plot, the lines are the empirical cumulative function of the trait values: blue solid lines refer to the trait values that lead to metric values below the metric-specific threshold; the red dashed lines refer to trait values that lead to metric values above the threshold. Shaded areas identify ranges of trait typically observed in mid-to-high latitudes. In the insets, the box and whisker plots compare and summarize the trait values leading to metric values above (left box; A–) and below (right box; B-) the threshold. The boxes extend from the first to the third quartile, the end of the whiskers are computed as 1.5 × IQR (Interquantile Range). The median values are indicated by the horizontal red lines.
