Quantifying in situ phenotypic variability in the hydraulic properties of four tree species across their distribution range in Europe

Many studies have reported that hydraulic properties vary considerably between tree species, but little is known about their intraspecific variation and, therefore, their capacity to adapt to a warmer and drier climate. Here, we quantify phenotypic divergence and clinal variation for embolism resistance, hydraulic conductivity and branch growth, in four tree species, two angiosperms (Betula pendula, Populus tremula) and two conifers (Picea abies, Pinus sylvestris), across their latitudinal distribution in Europe. Growth and hydraulic efficiency varied widely within species and between populations. The variability of embolism resistance was in general weaker than that of growth and hydraulic efficiency, and very low for all species but Populus tremula. In addition, no and weak support for a safety vs. efficiency trade-off was observed for the angiosperm and conifer species, respectively. The limited variability of embolism resistance observed here for all species except Populus tremula, suggests that forest populations will unlikely be able to adapt hydraulically to drier conditions through the evolution of embolism resistance.


Introduction
Massive forest mortality events due to drought stress and rising temperatures have been observed at the global and regional scales [1][2][3][4][5]. Considering that climate change models predict further increases in mean temperature and in the frequency and severity of extreme drought events [6], more negative impacts on tree survival are expected [7]. In this context, a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 assessments of the properties associated with drought resistance in trees, and of the capacity of species to deal with environmental changes, may help us to anticipate the impact of climate change on forest tree species.
Hydraulic failure due to xylem embolism is now considered one of the main causes of drought-induced tree mortality [8][9][10]. When soil water potential drops due to water shortage, the tension of the xylem water column increases, promoting the formation of embolisms that reduce the hydraulic functioning of the plant [11,12]. In cases of prolonged drought, soil water potential continues to fall, triggering the spread of embolisms throughout the xylem conduit network, leading to the hydraulic dysfunction of the plant vascular system and, finally, to lethal damage to the plant [13,14,15]. Therefore, determining the resistance to embolism of the species is crucial for evaluating the consequences that the expected increase in drought event frequency can have on a given population, forest or biome. P 50 is the xylem pressure at which 50% of conductivity is lost due to embolism formation, and it is widely used to assess plant hydraulic safety to embolism. Xylem-specific hydraulic conductivity (K S ), i.e. the rate of water transport through a given area of sapwood per unit pressure difference and per unit length, is commonly used to assess hydraulic efficiency [12]. Across species, literature shows a weak correlation between hydraulic safety and hydraulic efficiency, but the absence of species displaying both high hydraulic efficiency and safety suggests a possible safety-efficiency trade-off [16]. In conifers, P 50 and K s are only weakly correlated, as embolism resistance is driven mostly by the torus-aperture overlap in pit pairs [17][18][19], whereas xylem hydraulic efficiency is not influenced by this pit trait. By contrast, in angiosperms, both P 50 and K s are associated with pit membrane structure [20][21][22] and thickness [23], as well as with the perforation structure [24,25].
Differences in resistance to embolism, i.e. in P 50 , across species have been widely reported [17,19,26,27,28]. However, less attention has been paid to within-species phenotypic variation in this hydraulic property. Phenotypic variability results from a combination of genetic variation (differences in genotype among different individuals within the population and between populations) and phenotypic plasticity (genotype property to render different phenotypes in different environments [29]), and defines the capacity of populations to succeed under changing environmental conditions [30,31]. Low levels of phenotypic variability across large spatial scales may indicate a low potential of species to adapt to ongoing climate change. Contrary to other plant functional properties (see for instance [32] for leaf phenology, [33] for leaf functional traits), previous works on hydraulic properties show that phenotypic differences within species are by far lower than those found across species [34][35][36][37][38], and these differences are even smaller in gymnosperms than in angiosperms [39]. For instance, Lamy et al. 2011Lamy et al. , 2014 found neither phenotypic variability in situ and nor genetic differentiation between maritime pine populations and suggested uniform selection rather than genetic drift, for P 50 . However, whether the phenotypic variation of hydraulic properties varies across species distribution ranges remains largely unexplored. Furthermore, studies assessing the extent to which phenotypic variability in hydraulic properties is lower than that of other key species traits over large scales are also lacking.
The main aim of this study was to evaluate phenotypic variability in the functional hydraulic safety and efficiency (P 50 and K S ) of four European tree species (two conifers and two angiosperms) along a latitudinal gradient covering most of their distribution range. We assessed the capacity of the species to adapt to changing environmental conditions by exploring the links between hydraulic properties and latitude and climate. We also evaluated the phenotypic variability of branch growth in trees from the same populations to assess the extent to which hydraulic properties were conserved relative to other key traits. Finally, we assessed the safetyefficiency trade-off at intraspecific level. We hypothesized 1) a weaker phenotypic variation for hydraulic safety than for hydraulic efficiency and branch growth, given the highly conserved evolutionary nature of P 50 [37,38]; 2) a phenotypic cline -a gradual change of a phenotypic character in a species over a geographical area-in both hydraulic safety and efficiency and 3) a weak safety-efficiency trade-off within species. This study provides for the first time a multispecies assessment of inter and intra-specific phenotypic variability in functional hydraulic properties along a large latitudinal gradient. Our results will help to characterize the adaptive capacities of European forests, which will have to face drier and warmer climatic conditions in the future.

Study species and populations
We focused on four widely distributed European species, with different water-transport structures, from diffuse porous with scalariform perforation plates (Betula pendula Roth) or with simple perforation plates (Populus tremula L.) to softwood (two tracheid-bearing species, Picea abies (L.) Karst and Pinus sylvestris L.). For each species, four to six populations were selected across their distribution range (see the distribution range of the species and the location of the populations in S1 Fig). The mean annual temperature and total annual rainfall across selected populations ranged from -1.8 to 9.5˚C and 538 to 1739 mm, respectively (S1 Table). We also selected two different sites a few kilometres apart, for each population.

Climatic data
Data for mean annual temperature (MAT) and total annual precipitation (MAP) were obtained from WorldClim original 30-s data (http://www.worldclim.org/bioclim) [40] downscaled to 100-m resolution based on a high-resolution digital elevation model (DEM) and moving window regression technique [41] for all but the Italian (IT) and Swiss populations (SW-LOE and SW-PFY). The MAP data for the IT (Italy), SW-LOE (Switzerland-Loetschental) and SW-PFY (Switzerland-Pfynwald) populations were obtained from nearby weather stations at San Vito di Cadore (Centre for Alpine Environment Studies) and Sierre (www. meteoswiss.ch), respectively, due to considerable variations in topography. The aridity index (AI) was calculated as MAP/PET (total annual precipitation/annual potential evapotranspiration). PET was extracted from the Global Aridity and PET Database (http://www.cgiar-csi. org). We averaged the mean temperatures (T_Sum) or aridity indices (AI_Sum) of June, July and August to obtain mean values for the summer (see S1 Table for the climatic conditions of the populations studied).

Xylem vulnerability to embolism
We collected branches from five to 11 healthy mature trees per population in the early morning during the wet season (spring 2015). One or two branches with three to five functional rings were sampled at mid-crown and south oriented. Samples had a standard length of 45 cm. Transpiration losses were prevented by removing the leaves or needles immediately after sampling and wrapping the branches in moist paper to keep them humid and cool (3˚C) until the measurement of embolism resistance (within three weeks of sampling). The bark was removed from conifer branches to prevent resin to fill the cavitron reservoirs (see below, [17]), and all branches were recut with a razor blade, under water, to a standard length of 0.27 m. For each angiosperm species, 10 samples per species were used to test the open vessel artefact [42] by injecting air at 2 bars at one end and no open vessels were detected for any of them.
Vulnerability to drought-induced embolism was determined at the Caviplace (University of Bordeaux, Talence, France; http://sylvain-delzon.com/caviplace) and INRA-Clermont-Ferrand facilities, with the Cavitron technique [43,44]. Samples were infiltrated with a reference ionic solution of 10 mm 25 KCl and 1 mm CaCl 2 in deionized ultrapure water. Centrifugal force was used to generate negative pressure into the xylem and induce cavitation. This method allows to measure xylem conductance under negative pressure using the custom software Cavisoft 4.0 (Univ. Bordeaux, Pessac, France). Initially, the maximum conductance of stem (K max , in m 2 MPa -1 s -1 ) was calculated under low xylem pressures. The percentage loss of conductance (PLC) of the stems was calculated at different xylem pressures (P i ) from -0.8 to -5 MPa with the following equation: We obtained one vulnerability curve per tree by measuring one or two of the collected branches. These vulnerability curves show the percentage loss of xylem conductance as a function of xylem pressure [17]. For each branch, the relationship between PLC and xylem water pressure was fitted with the following sigmoidal equation [45]: where P 50 (MPa) is the xylem pressure inducing a 50% loss of conductivity and S (% MPa -1 ) is the slope of the vulnerability curve at the inflection point. All sigmoidal functions were significant and fitted with the NLIN procedure in SAS (version 9.4 SAS Institute, Cary, NC, USA). The xylem-specific hydraulic conductivity (K s , kg m 1 MPa -1 s -1 ) was calculated by dividing the hydraulic conductivity measured at low speed by the sapwood area of the sample.

Branch growth measurements
We also collected one branch per tree from three to five trees per population and site. We selected straight branches and did not keep any sample with reaction wood for our measurements. The allometric relationship between branch radius (mm) and xylem age (number of years) was used as a surrogate for tree growth, as radial branch growth and tree growth patterns are highly correlated [46]. The branch surface area and the number of tree rings were systematically measured at 70 cm from the branch apex.

Statistical analyses
We assessed the phenotypic variability of functional hydraulic properties (P 50 and K s ) and branch growth (branch radius/xylem age) in each species, by testing the effect of population and site with nested ANOVAs, in which the population and the site nested in population were considered factors. If statistically significant differences were observed, post-hoc Tukey tests were conducted for multiple comparisons between populations. Before running the ANOVAs, we checked that the data satisfied the assumptions of normality and homoscedasticity. As vessel size can rapidly increase with branch size during early years of tree growth, and then may have a potential effect on hydraulic conductivity, we tested any potential correlation between K s and branch diameter. We also calculated the inter-population and intraspecific coefficients of variation (% CV inter and CV sp , respectively). For each species, Spearman´s or Pearson´s correlation coefficients (depending on the linearity condition) were calculated between the averaged by site P 50 , K s and branch growth and the latitude and the five climatic variables mentioned above. Finally, we also checked for intraspecific safety-efficiency trade-offs, with Spearman´s or Pearson´s correlation tests. Statistical analyses were performed with the R project for statistical computing (R Development Core Team, 2016) [47].

Phenotypic variability across species distribution ranges
Xylem vulnerability curves followed a sigmoid function in all species (Fig 1, S2 Fig), showing the lack of an open vessel artefact and the accuracy of the results obtained here. Betula pendula and Picea abies showed, respectively, the lowest and highest resistance to embolism of the four species evaluated. The mean P 50 ± SE (MPa) was -1.78 ± 0.02 for Betula pendula, -2.45 ± 0.08 for Populus tremula, -3.16 ± 0.03 for Pinus sylvestris, and -3.58 ± 0.02 for Picea abies (Fig 1).
Differences in P 50 between populations were observed for all species, whereas the site (nested in population) had an effect on P 50 in all species but Picea abies ( Table 1, Fig 2). The CV in P 50 was low for all species other than Populus tremula. The variability in P 50 of Betula pendula, Picea abies and Pinus sylvestris ranged from 4.15 (CV inter of Picea abies) to 10.23% (CV sp of Pinus sylvestris), whereas that of Populus tremula ranged from 24.82 (CV inter ) to 25.07% (CV sp ) ( Table 2). The high variability observed for Populus tremula was mostly due to the population of Finland-Värriö (FI-VA), which had the least negative P 50 values of any of the populations studied (Fig 2). Variability levels were much higher for K s than for P 50 ( Table 2, Fig 2). K s differed significantly between populations, for all species other than Populus tremula ( Table 2, Fig 2). For this    Fig 2). The phenotypic variability of branch growth was greater than that of P 50 (Table 2, Fig 2). Furthermore, the phenotypic variability of branch growth was greater than that of K s in most cases ( Table 2, Fig 2). Pinus sylvestris had the largest CV sp and CV inter in branch growth (85.81 and 76.67%, respectively), whereas these two coefficients were the lowest in Populus tremula (43.29 and 37.66%, respectively) ( Table 2).

Phenotypic clines with climate and latitudinal gradients
Populus tremula presented strong significant clines in P 50 , as five out of the six climatic variables studied here showed significant correlations with P 50 (Table 3, Fig 3). P 50 values for this species were positively correlated with latitude and aridity index (AI and AI_Sum), but negatively correlated with MAT and T_sum (Table 3, Fig 3). There was also a statistically significant negative correlation between P 50 and T_Sum in Betula pendula (Table 3). By contrast, no significant clines in P 50 were observed for conifers (Table 3, Fig 3).
K s was less strongly related to climate than P 50 . The K s /climate correlation was statistically significant only for Betula pendula, with lower K s values at sites with higher MAT values ( Table 3, Fig 4).
We found steeper clines for branch growth than for hydraulic properties, with all species showing at least one statistically significant correlation between branch growth and latitude/ climate variables (Table 3, Fig 5). Latitude was correlated with branch growth only in Populus tremula, for which the lowest branch growth values were obtained for the northernmost population (Table 3, Fig 5). In general, when statistically significant, branch growth was positively correlated with MAT and T_sum, and negatively correlated with AI, AI_sum and MAP (Table 3, Fig 5).

Safety-efficiency trade-off
At the intraspecific level, we found statistically significant but weak positive correlations between P 50 and K s for the conifers studied, with the most vulnerable individuals having the Table 2 50 and K s for either of the angiosperms studied (Table 4).

Discussion
We assessed the phenotypic variability of hydraulic safety and efficiency traits (P 50 and K s , respectively) and branch growth in four tree species across a long latitudinal gradient covering most of their distribution range in Europe. P 50 displayed lower phenotypic variability than K s and branch growth, consistent with our initial hypothesis. The low variability of P 50 across populations has been related to uniform evolutionary selection or canalization [37,38]. Indeed, these studies provided evidence of natural selection acting on this trait. This uniform selection reproduces trait conservatism and eventually leads to stasis [37]. In contrast, it has been suggested that K s variability is related to the interaction between genotype and environment [35]. K s may also vary significantly with sampling position along the branch axis [48,49], although we tried to overcome this limitation through the use of systematic sample preparation procedures. We also expected branch growth to be more variable than P 50 , because branch growth is strongly influenced by multiple interacting factors, such as the availability of nutrients, light, water and temperature [46,50,51,52,53] and biotic interactions [54]. The limited embolism resistance variability observed here, in all species other than Populus tremula, suggests that forest populations of the studied species will potentially find difficulties to cope with a warmer and drier conditions by increasing their embolism resistance. However, considering the differences in climatic ranges between the studied species, we have to be cautious when interpreting these patterns. Further studies investigating larger precipitation gradients and/or marginal Table 3 populations are needed. Indeed, a recent study showed that marginal populations of beech significantly differed in embolism resistance [55] while core populations exhibited similar P 50 values [36]. Our results also show that within species phenotypic variability in K s and growth are large, and in general larger than that of P 50 , suggesting that intra-population variability should not be neglected in further studies at local scales. Phenotypic variability in hydraulic traits Across species, higher resistance to embolism (i.e. more negative P 50 values) have been related to drier environments [26,27]. However, at the intraspecific level, no statistically significant correlations between P 50 and climate have been observed for herbaceous plants [56], angiosperm trees [57,58] and conifers [34,38,59,60]. Our results are consistent with these findings, as only one species, Populus tremula, presented a significant cline in P 50 . However, the potential effect of collinearity between climate variables could not be investigated here. Previous studies have reported higher [59,61], similar [34,60] or lower K s values [35,57] at dry sites than at mesic sites. However, we found no support for a K s /climate cline. In contrast to P 50 and K s , all species showed significant correlations between climate and growth. Populations growing at high latitudes and in cold temperatures had the lowest levels of branch growth, probably due to the shorter growing season. Phenotypic variability in hydraulic traits Ideally, plants should be able to maintain both efficient conductivity and the safety of the hydraulic system. However, this is not always the case in natural conditions, and little or no support for a safety-efficiency trade-off has been obtained across species [16]. A few studies have evaluated this trade-off at intraspecific level, and found either no support for the existence of a trade-off [34,62] or an association of greater conductivities with lower embolism resistance [63]. The lack of correlation between P 50 and K s in angiosperms, and the weak correlation found here for conifers provide insufficient support to conclude that there is a safetyefficiency trade-off. Current knowledge of the anatomical basis of P 50 and K s in conifers is also consistent with the absence of support of such a trade-off, as P 50 is determined principally by the torus-aperture overlap in this clade, whereas K s is not related to this anatomical trait [17][18][19] but rather to the vessel lumen area.

Betula pendula
Researchers have recently developed an interest in the phenotypic variability of hydraulic properties [34,35,37,39], due to its possible contribution to community assemblages. However, studies of phenotypic variability can also reveal the potential of the species to adapt to the new environmental conditions imposed by ongoing climate change. We found statistically significant differences in hydraulic safety traits between populations in all the species studied here, and a lack of climate cline in all species other than Populus tremula. Overall, this species displayed the greatest variability of P 50 within species and between populations, suggesting a potentially greater ability to adapt to environmental changes. This higher variability in safety traits might be due to the fact that Populus tremula frequently hybridizes with Populus alba in Europe [64]. Yet, the two species significantly differ in numerous phenotypic and ecological properties [65]. In contrast, the low intraspecific variability and lack of climate clines for hydraulic properties in conifers suggest strong genetic constraints, with a much smaller potential to evolve greater embolism resistance in the xylem to cope with the predicted drier conditions. However, hydraulic adjustments can also occur through changes in leaf area: sapwood area ratio [34,66], and a decrease in transpiring leaf area relative to xylem conductive area could hence maintain plant water balance under drought conditions.

Conclusions
Quantification of the adaptive capacity of populations and species is important for the prediction of natural adaptation to climate change, especially in the long term. Adaptation requires the presence of genetic variation among the individuals of populations upon which natural selection can act. The phenotypic variability in embolism resistance found here was weak and much smaller than that for xylem conductivity and branch growth. In addition, no relationship was found between embolism resistance and climatic variables, except for Populus tremula. The species studied (except Populus tremula) would therefore be unlikely to be able to adapt hydraulically to drier climatic conditions through the evolution of embolism resistance. Our results provide little support to the existence of a hydraulic safety-efficiency trade-off at the species level.