The plant surface temperatures can be related to drought stress and plant water status by several methods (reviewed in Maes and Steppe, 2012). To account for the changing environmental conditions, the surface temperatures must be normalized by reference temperatures (see Reference surfaces and plant monitoring platform). We used two thermal indices; the crop water stress index (CWSI) and stomatal conductance index (Ig).

Jones (1999) proposed the CWSI as a modification of Idso et al.'s (1981) formulation. The CWSI is known to mirror stomatal conductance and the leaf and stem water potentials. It normalizes the leaf surface temperatures by the surface temperatures of wet (T_wet) and dry (T_dry) references, where T_wet represents a fully transpiring leaf and T_dry a non-transpiring leaf. The CWSI is calculated as (1)CWSI=Tcanopy - TwetTdry - Twet

I_g employs the same variables as CWSI but is linearly related to the stomatal conductance (Jones, 1999). Thus, Ig is a linear function of the stomatal opening: (2)Ig=Tdry - TcanopyTcanopy - Twet

Higher CWSI and lower Ig values indicate higher surface/tissue temperatures and thus stomatal closing.

As suggested in Meron et al. (2003) and Möller et al. (2006), we calculated the abovementioned thermal indices using artificial reference surfaces (ARSs), which are included in each picture (see Meteorological data and Figure 2). To mimic the maximum transpiring surface, we wrapped a white cotton fabric around a styrofoam board floating in a water-filled plastic box (wet reference). A non-transpiring leaf was represented by an opal white laminated fiberboard (dry reference) mounted on the plastic box (Figure 2). Jones (1999) used wet and dry leaves as the reference surfaces, but here we chose the ARSs because the thin needles of conifers (unlike leaves) easily dry out under the high greenhouse temperatures (Figure 3A). Additionally, the temperature information in the pixels of the reference needles might become mixed with that of non-reference needles in the background. Unfortunately, this approach might be sensitive to changing environmental conditions because the short-wave absorptances and heat capacities differ between needles and reference targets. For these reasons, environmental variables were incorporated as control covariates in the statistical models (see Statistical analysis). In addition to the wet and dry ARSs, two black painted electric heating plates with a mean temperature of 40°C were horizontally attached to a support frame. These created a strong contrast to the plant tissue and totally masked the pot and soil in the images. A small gap in the middle edge of the plates prevented squeezing of the trunk. The handling time of taking an individual tree to the monitoring platform, mounting it into the platform and capturing up to three thermal images was approximately 3 min.

Thermal images were acquired by a thermal infrared camera (VarioCam hr inspect 780, Infratec, Dresden Germany) with a resolution of 1280 × 960 infrared pixels using the optomechanical resolution enhancement. The thermal resolution within the images was below 0.08 K at 30°C and the absolute measurement accuracy was ± 1.5 K. The emission coefficient was set to 1 during the image acquisition and was later corrected using the emission coefficients calculated by the quotients of the contact thermocouple temperatures and the thermal image-derived temperatures of the reference surfaces and pine seedlings (performed in a dark chamber under temperature and relative humidity control and in the greenhouse, respectively). The resulting emission coefficients of the wet reference, dry reference and pines were 0.95, 0.93, and 0.90, respectively. The calculated emissivity can be considered as the apparent emissivity since we did not measure background temperature. The errors caused by this omission should be negligible inside the greenhouse, which was constantly shaded for the thermal imaging, largely corresponding to cloudy conditions (see Meteorological data, Maes and Steppe, 2012). Thermal images were taken approximately twice a week from July 5th to September 4th, leading to 17 measuring days. Ten of these days constituted the stress period with clear PSWD effects, and 4 days followed the summer drought treatment (see Figure 1B for exact dates). During the image acquisition (between 11 a.m. and 3 p.m.), the greenhouse was shaded to reduce the possible influence of variable solar radiation (e.g., due to scattered clouds). The camera was vertically mounted at 2.5 m above the monitoring platform attached to the scaffolding of the greenhouse. On every measuring day, each of the 48 pines was photographed 2 or 3 times (rarely once due to technical problems), together with the reference surfaces.

In the subsequent image analysis, the plant tissue was separated from the (heated) background, and interferences at the needle edges (where single pixels were mixtures of plant and background temperatures) were additionally removed. To automate the image processing, a script was written in Fiji (Schindelin et al., 2012). Images were processed by the following steps. First, the image was sharpened using the command “Unsharp Mask” with a radius of 2 and a mask weight of 0.9. From this sharpened image, two masks were created, one to remove the background (using auto thresholding based on the intermodes algorithm), the other to remove the edges (using the “Find Edges” command). The median temperature of the plants' canopy and the dry and wet references in each image were calculated from the remaining pixels. The respective mean daily temperatures were then determined from the multiple (1–3) images acquired on each measuring day. Finally, the thermal indices were calculated from the mean temperatures per measuring day and used in subsequent statistical analysis.

During the thermal image acquisition, an air temperature (T) and relative humidity (RH) sensor (HOBO U23 Pro v2, Hobo®, Onset Computer Corporation, Bourne, MA) was placed next to the plant monitoring platform. Data were recorded in 1-min intervals. In the greenhouse, the air temperature, relative humidity and solar radiation were measured at 10-min intervals throughout the whole study period by a meteorological weather station (Davis Vantage Pro2 Plus™, Davis Instruments, Hayward, CA). We matched the meteorological and thermal image data with their nearest temporal counterparts. The vapor pressure deficit (VPD), defined as the difference between saturation vapor pressure (e_s) and actual vapor pressure (e_a), was calculated after Allen et al. (1998) with T and RH as the input variables.

The daily mean values of the air temperature, relative humidity, solar radiation, and vapor pressure deficit, collected at 10-min intervals between 11 a.m. and 3 p.m. (the time window of the thermal image acquisition on measuring days), varied during the study period (July 5th to September 4th), within the ranges 16.6–36.6°C, 26.4–85.0%, 33.0–421.7 W/m², and 0.3–4.5 kPa, respectively (Figure 3). The minimum and maximum air temperature, relative humidity, solar radiation and vapor pressure deficit during the image-acquisition time over the 17 measuring days were 20.4 and 37.1°C, 29 and 66%, 24 W/m² and 234 W/m², and 1.0 kPa and 4.5 kPa, respectively. Thus, the conditions during measuring times well represented the indoor conditions at noon over the course of the study period except for the solar radiation, which was mostly determined by the additional shading in the greenhouse (generally throughout June and July, and on the measuring days after mid-August).

We analyzed the effect of PSWD on thermal indices over the whole study period (from July 5th to September 4th). The analysis was performed by linear mixed-effects models (nlme; Pinheiro et al., 2016) implemented in R version 3.2.2 (R Core Team, 2015). Full models were constructed by adding the covariates provenance as factorial dummy variable, heights, and diameters of seedlings at the beginning of the experimental period, air temperature, relative humidity, vapor pressure deficit, and solar radiation and the two-way interactions of PSWD with provenance, air temperature, relative humidity, and vapor pressure deficit. Because the PSWD varies nonlinearly with the thermal indices, it was added as a linear and quadratic term to the models. All covariates besides PSWD were centered on their means (by subtracting their respective means from the discrete variable values) for easier interpretation of interaction effects. As time and PSWD were strongly correlated, they cannot be included simultaneously due to collinearity; hence, no time variable was included in the PSWD-based models. To account for the repeated measurements of individuals during the experiment, we included the individual trees as random variables.

In a first step of model selection we either chose air temperature plus relative humidity or vapor pressure deficit as covariates based on the AIC (Akaike Information Criterion) of the respective models, since these variables were strongly correlated. The full model is mathematically expressed as (3)Indexi,j=β0 + b0,i + γ1Provenancei + β1PSWDi,j                      + β2PSWDi,j2 + β3Heighti + β4Diameteri                      + β5VPDi,j + β6Radiationi,j + γ2Provenancei                      * PSWDi,j + γ3Provenancei * PSWDi,j2 + β7VPDi,j                      * PSWDi,j + ϵi,j with random intercepts b0,i~N(0,σ02) and errors ϵi,j~N(0,σϵ2). Index_{i, j} is the thermal index of tree i (= 1 …48) at measurement j (= 1 …17). The meteorological variables (in this case, the VPD and radiation) also vary from tree to tree, because they were measured sequentially on the measurement day, and altered throughout the course of the day. As the six provenances were modeled using dummy variables, each of γ₁, γ₂, andγ₃ is a five-dimensional coefficient vector of the dummy regressors.

To simplify the full models, we evaluated the importance of the explanatory variables/interactions using the drop1 function (stats; R Core Team, 2015). Any terms that did not improve the models' explanatory power were excluded (Table 2). The R² of the final models was computed by the r.squaredGLMM function (MuMIn; Barton, 2015, Table 3). To avoid heteroscedacity and non-normal distribution of the residuals, we examined the diagnostic plots and applied variance function structure classes. The Ig was square-root-transformed to meet these criteria. Provenances were compared by a pairwise post-hoc test using the glht function (multcomp, Hothorn et al., 2008) comparing contrasts with the Tukey's range test. By this test, we also compared the provenances under additional PSWD scenarios (0%, 50%, and 100% PSWD) after centering the PSWD values on these thresholds and refitting the models. To check the different stress behaviors of provenances in relation to PSWD, we tested the coefficients of the linear and quadratic terms of the PSWD–provenance interaction with the glht function.

Differences in thermal indices between the treatments (water supplies) and among the provenances, and in their corresponding ability of the pines to recover from the water stress, were separately analyzed over the stress period (measuring days 4–13, July 17th until August 21st) and the recovery period (measuring days 14–17, August 24th until September 4th). Here we fitted and simplified the linear mixed models as described above. The factorial dummy variable treatment (control/drought) was used rather than the PSWD, and the two-way interactions of treatment with height and diameter were added to the initial full model. As the PSWD changes over time in both treatments and study periods, we added a time variable (number of days since the observations started) as a control covariate with a linear and a quadratic term. Based on the AIC, we selected the VPD and air temperature/relative humidity as covariates. The full model is mathematically expressed as (4)Indexi,j=β0 + b0,i + γ1Provenancei + β1Treatmenti + β2tj                 + β3tj2 + β4Heighti + β5Diameteri                + β6Temperaturei,j + β7RelativeHumidityi,j                + β8Radiationi,j + γ2Provenancei * Treatmenti                + β9Treatmenti * tj + β10Treatmenti * tj2                + β11Treatmenti * Heighti+ β11Treatmenti                * Diameteri + β12Treatmenti * Temperaturei,j                + β13Treatmenti * RelativeHumidityi,j + ϵi,j where Treatment_i is 0 if tree i received the dry treatment, and 1 for the control group. t_j is the number of days after the start of the experiment (July 5th) at measurements j (= 4 …13) for the stress period and j (= 14 …17) for the recovery period. Employing the glht function, we again tested the differences among provenances and between treatments in pairwise post-hoc tests, comparing contrasts with the Tukey's range test.

Additionally, we analyzed the response magnitudes of the provenances between treatments during the stress period. The mean index values were calculated for each individual. The response magnitude was calculated as the pairwise difference in particular indices between each individual of the drought group and all other individuals of the control group. Differences in response magnitude, tree height and diameter were analyzed by the Kruskal–Wallis test (stats; R Core Team, 2015) and by Dunn's test for multiple comparisons (FSA; Ogle, 2015).

All p-values for multiple comparisons were corrected by the false discovery rate (FDR) method.

The two thermal indices were significantly related to PSWD, provenances and meteorological control parameters. Excluding some meteorological control covariates (temperature, relative humidity) and plant biomass covariates (diameter, height) from the full models, the final models resulted in better fits in terms of lower AIC-values and explained a proportion of variance of 0.64 and 0.62 for CWSI and Ig, respectively (Table 2).

To interpret the signs of the model variables (Table 2), we must remember that CWSI and Ig increase and decrease at higher stress levels, respectively. In both thermal indices, higher stress levels were linked to higher PSWD, indicating lower soil-water availability (Table 2; Figure 4). However, provenance and the its interaction with PSWD were also included in the models of the thermal indices (Table 2). The Bulgarian and especially the Spanish provenance showed a pronounced stress minimum (PSWD = 15–20%; see Figure 4). Although responses of provenances to PSWD were varying, significant differences were detected only in the Spanish provenance. Pines in this provenance were much more sensitive to increasing water deficit than all provenances other than the Bulgarian one (Table 3; Figure 4). Increasing VPD and radiation also increased the stress levels (either by natural response of the stomata and leaf surfaces or because the ARS-derived indices were sensitive to the changing environmental conditions); Figure 4 displays the responses of the different provenances to PSWD, for the mean VPD and radiation conditions. Note that all dimensional variables (seedling height and diameter) were excluded from the final model because they added no further explanatory power aside from provenance.

Under well-watered conditions (PSWD = 0%) the provenances from Spain and Bulgaria were less stressed than those from Poland and Germany (p < 0.05; Supplementary Table 1). Marginal differences were observed between pines from France and Germany and between pines from Italy and Spain (p < 0.1; data not shown). At medium soil moisture (PSWD = 50%) the Bulgarian and Spanish seedlings were less stressed than the German and Polish seedlings and the Spanish provenance performed better than the Italian provenance (p < 0.05). Additionally, the CWSI was marginally different between the French and Spanish seedlings and between Italian and Bulgarian provenances (p < 0.1; data not shown). Under severe drought conditions (100% PSWD), only the two most extreme provenances (France and Spain) showed significant differences in both thermal indices (p < 0.01). The Ig significantly differed between pines from Poland and France (p < 0.05; Supplementary Table 1).

The immediate drought stress during the stress period and the conditions during the recovery period, were assessed for both treatments (defined by their water supply i.e., the amount of irrigation water) by the thermal indices. For both CWSI and Ig, the explanatory power of the corresponding linear mixed models was lower during the stress period (both indices = 0.63) than during the recovery period (0.74 and 0.75, respectively, Table 4).

During the stress and recovery periods, the respective final models of CWSI and Ig included provenance, time and treatment. As expected, stress was increased by the drought treatment and continued to increase throughout the stress period, then declined over the recovery period. In contrast, the sets of plant specific (dimension, provenance) and meteorological covariates explaining the variation in thermal indices differed between stress and recovery periods (Table 4).

Diameter was negatively associated with stress during the stress period, but was unimportant in the models of the recovery phase (Table 4; Figures 5B,C). The influence of seedlings' height on the thermal indices depended on the water supply treatment (Table 4; Figures 5B,C). Specifically, the stress levels increased with height in the control treatment but were independent of height in the drought treatment. Equally, height was unimportant in the recovery period models.

During the stress period, the effect of provenance on the thermal indices depended on the treatment (Table 4; Figure 5), whereas during the recovery period, the treatment and provenance influenced the thermal indices without interacting with each other (Table 4; Figure 6). In each provenance, the seedlings were significantly more stressed in the drought than in the control group during the stress period (p < 0.01, Supplementary Table 2). The exception was the Italian provenance, where high stress levels were already observed in the control group (Figure 5). Although the treatment differences under stress decreased by ~50% during the recovery period (Table 4), the overall differences between treatments persisted in the recovery period (Figures 6B,D; Supplementary Table 3). During the stress period, the provenance differences were not significant in the highly stressed drought treatment group. However, in the control group, whose PSWD levels imposed moderate stress during this period, the CWSI levels indicated more stress in the Italian provenance than in provenances from France, Spain and Bulgaria. However, the Ig did not capture these differences at a significance level of 0.05 (Figures 5A,D; Supplementary Table 2).

During the recovery period, provenances showed equal stress patterns in both treatment groups, because the provenance and treatment interaction was not significant in the statistical model (Table 4). According to both thermal indices, the provenances from Spain and Bulgaria were less stressed than those from Poland and Germany, and the French provenance was less stressed than the German provenance (Figures 6A,C; Supplementary Table 3). Furthermore, the CWSI showed higher stress levels in the Italian seedlings in than the Spanish and Bulgarian ones. Although the seedlings recovered, the significant difference between the control and experimental groups persisted into the recovery period, indicating some post-drought stress (Table 4; Figures 6B,D).

The treatment effect (water supply) during the stress period clearly separated the individuals of almost every provenance. The exception was the Italian provenance, in which the negative group differences indicated overlap of the stress levels of both treatment groups (Supplementary Figure 3). The differences between treatments exhibited similar patterns for both thermal indices. The Spanish and Italian provenances demonstrated the strongest and weakest response to water supply, respectively. During the stress period, the response to the water supply treatment was stronger in the Spanish pines than in the Polish, Bulgarian, and Italian pines, as shown by both indices (Supplementary Table 4). The CWSI and Ig revealed further differences in the water supply response; specifically, the Italian seedlings demonstrated weaker CWSI response than the French, Polish, and German seedlings, and the Ig response to water supply was higher in Spanish than in German pines (Supplementary Table 4).

When performing thermal imaging under non-100% controlled conditions, several limitations must be considered (Prashar and Jones, 2014). Variations in environmental conditions affect the leaf temperature, so reference surfaces are required. Here we replaced wet and dry Scots pine needles with ARSs (see Materials and Methods) instead of using wet and dry Scots pine needles which would have had similar spectral and physical properties to the target leaves (Jones, 1999). Nevertheless, the ARS performance should be adequate for the following reasons: (1) the emissivity of each target surface was determined; (2) the greenhouse was shaded before and during image acquisition, meaning that all thermal images were acquired under low levels of incoming solar radiation; (3) on most measurement days, the T, RH, and VPD varied within a narrow range, suggesting that the errors introduced by differing heat capacities of the targets are also small; (4) when comparing indices across provenances and treatments, these errors can be partially compensated by incorporating meteorological covariates in the statistical models.

Additionally, the temperature within a thermal image might vary with angle of view, illumination, and distance to targets (Prashar and Jones, 2014). In our experimental design, these parameters were maintained largely constant by a plant monitoring platform. All images were horizontally captured above the seedlings during a short time interval, ensuring a fixed distance to the seedlings. The mixed pixels on the plant edges were removed during the image processing. As all seedlings were morphologically similar, these geometric and lighting factors can also be neglected.

The differences in stress sensitivity and recovery of the six provenances were studied after adjusting for PSWD or (more simply) for treatment (water supply). During the stress period, the increased PSWD indicated mild drought stress even in the control group (Figure 1B). Therefore, the provenance performance was studied under two conditions of water stress; extreme (drought treatment) and moderate (control treatment in the stress period).

The Spanish and Bulgarian seedlings were more distressed by over-saturated than under-saturated water conditions (Figure 4). Excessive soil water may decrease the oxygen availability for roots. The lower root respiration reduces the root function and causes leaf dehydration (Vartapetian and Jackson, 1997). Consistent with these findings, the ecophysiological leaf traits of Scots pine (net assimilation, stomatal conductance and transpiration) are lower in short-term water-logged soils than in soils hydrated to field capacity (Repo et al., 2016). The germination rates, shoot growths and root growths of Scots pines of different provenances also respond differently to waterlogging (Mukassabi et al., 2012).

Important triggers of the drought stress response besides PSWD were treatment, provenance and the provenance–PSWD and provenance–treatment interactions. However, because of the limited sample size, significant differences in stress sensitivity between the provenances were rarely identified. Both lines of evidence (PSWD and water supply) imply that pine species' response to strong drought stress (either modeled for 100% PSWD or for the drought treatment during the stress period, when the soil moisture approximated the wilting point) is comparable among provenances. At 100% PSWD, significant differences were observed only between the most extreme provenances (Spain and France, and additionally between France and Poland in the Ig analysis). In contrast, for moderate stress (i.e., over the range of measured PSWD, for the given water supply in the control treatment, or as modeled for 50% PSWD) our study revealed significant differences. Pines in the Spanish provenance tended to be less sensitive to moderate drought stress than Italian pines. Furthermore, the Spanish and Bulgarian pines were significantly less stressed at 50% PSWD than the German and Polish ones, and the Italian provenance demonstrated higher response to mild stress under the control water supply than the French and Bulgarian provenances. Taeger et al. (2013a, 2015) investigated the stem diameter, stem length and respective relative growth rates in pines from these provenances, and reported similar differences among the provenances. The difference in stress levels between the control group (exposed to moderate stress) and the drought group (exposed to extreme drought stress) were smallest and largest in pines from Italy and Spain, respectively (see Supplementary Figure 3).

The Spanish, French, and Italian provenances experience a Mediterranean climate with minimal precipitation in summer; whereas the Bulgarian provenance (despite its similarly dry summer season) is continental (see Supplementary Figure 4 and Table 1). In contrast, the German and Polish provenances are clearly characterized as temperate-continental.

In moderate drought scenarios (50% PSWD or the control group during the stress period in), the pines from provenances with a summer precipitation minimum at their origin (Bulgaria, France, Spain) were more stress-resilient than pines from provenances with rainy summers (Germany, Poland). Thus, the temperate continental provenances exhibited the highest stress under moderate drought. The superior resilience of the Mediterranean type group (the Spanish and Bulgarian provenances) lacks an obvious explanation. Mediterranean Scots pines, which are considered to naturally adapt to drought events (Cregg and Zhang, 2001; Richter et al., 2012; Taeger et al., 2013a), may have adapted their stomatal control under our study conditions. This finding may also imply that southern (Mediterranean) provenances transpired less and used a smaller share of the offered water supply 3050 ml during 42 days of treatment. Tognetti et al. (1997) reported comparable results in drought-stressed Pinus halepensis, whose leaf conductance depends on the moisture content of its origin. Drought adaptation in P. sylvestris might be governed by lower investment in aboveground biomass and higher biomass allocation to roots (Cregg and Zhang, 2001; Taeger et al., 2015). However, identifying the mechanism of drought adaptation is beyond the scope of this study.

During the recovery period, pines from the Polish, German, and Italian provenances maintained significantly higher stress levels than their Spanish and Bulgarian counterparts. This trend might also reflect a differential stomatal control.

An important feature of our study was the overlap of stress levels in both treatments, despite clear PSWD differences between the treatments. This phenomenon, which appeared in some instances and was especially observed in pines from the Italian provenance, suggests a resource saving strategy as the stress levels were already high under mild drought conditions (control treatment). Such a strategy is supported by a previous mortality experiment, in which trees from the Italian provenance were least threatened by drought-induced mortality (Seidel and Menzel, 2016).

The seedling diameter did not significantly differ among the provenances, but the continental individuals were taller than the Spanish and French ones. Although the seedling dimensions were not included in the final linear mixed-effects models relating PSWD to stress levels, provenance-specific growth traits might contribute to the water responses of individual trees, as mentioned above. The exclusion of tree dimensions from the PSWD models indicates that dimensional traits do not purely drive the observed differences between provenances. This is also indicated by the fact that the water supply models included both provenance and tree dimensions.

In the water supply models, the seedling dimensions (diameter and height) significantly affected the stress response. In particular, the stress levels generally decreased with increasing tree diameter (Figures 5C,F). Trees can store a considerable amount of water in sapwood, whose volume is closely related to stem diameter (Meinzer et al., 2001). Indeed, a simulation study of Pinus sylvestris showed a strong relationship between stored water use and tree diameter. Stored water use can contribute up to 40% of the total daily transpiration (Verbeeck et al., 2007). Thus, large diameter provides a water buffer against drought stress. Under mild drought conditions, we detected a positive relationship between stress level and seedling height (Figures 5B,E). A link between reduced growth of aboveground biomass and increased drought adaptation has been suggested (Alía et al., 2001; Valladares et al., 2007). Thus, the increased stress levels in higher seedlings may reflect the higher water consumption of larger sized seedlings or the higher transpiration of larger needle area. In contrast, the seedlings in the drought group exhibited uniform stress levels, most likely because they consumed the irrigation water to an equal extent.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.