Research Methods Towards a statistically robust determination of minimum water potential and hydraulic risk in plants Jordi Martınez-Vilalta1,2 , Louis S. Santiago3 , Rafael Poyatos1,2 , Llorencß Badiella2 , Miquel de Caceres1,4 , Ismael Aranda5 , Sylvain Delzon6 , Alberto Vilagrosa7 and Maurizio Mencuccini1,8 1 CREAF, Bellaterra (Cerdanyola del Valles), Catalonia E08193, Spain; 2Universitat Autonoma de Barcelona, Bellaterra (Cerdanyola del Valles), Catalonia E08193, Spain; 3Department of Botany & Plant Sciences, University of California, 2150 Batchelor Hall, Riverside, CA 92521, USA; 4Joint Research Unit CTFC – AGROTECNIO, Solsona 25280, Spain; 5Centro de na Km 7.5, Madrid E-28040, Spain; 6INRAE, BIOGECO, Univ. Bordeaux, Investigacion Forestal, Instituto Nacional de Investigacion y Tecnologıa Agraria y Alimentaria, Carretera Coru~ Pessac 33615, France; 7CEAM Foundation, Joint Research Unit University of Alicante-CEAM, Dept Ecology, University of Alicante, Carr. de San Vicente del Raspeig, PO Box 99, Alicante 03080, Spain; 8ICREA, Pg. Lluıs Companys 23, Barcelona 08010, Spain Summary Author for correspondence:  Minimum water potential (Ψmin) is a key variable for characterizing dehydration tolerance Jordi Martınez-Vilalta and hydraulic safety margins (HSMs) in plants. Ψmin is usually estimated as the absolute mini- Email: Jordi.Martinez.Vilalta@uab.cat mum tissue Ψ experienced by a species, but this is problematic because sample extremes are affected by sample size and the underlying probability distribution. Received: 8 January 2021  We compare alternative approaches to estimate Ψmin and assess the corresponding uncer- Accepted: 14 June 2021 tainties and biases; propose statistically robust estimation methods based on extreme value theory (EVT); and assess the implications of our results for the characterization of hydraulic New Phytologist (2021) 232: 404–417 risk. doi: 10.1111/nph.17571  Our results show that current estimates of Ψmin and HSMs are biased, as they are strongly affected by sample size. Because sampling effort is generally higher for species living in dry Key words: drought resistance, drought environments, the differences in current Ψmin estimates between these species and those liv- tolerance, extreme values, hydraulic risk, ing under milder conditions are partly artefactual. When this bias is corrected using EVT meth- hydraulic safety margin, minimum water ods, resulting HSMs tend to increase substantially with resistance to embolism across species. potential, physiological limits.  Although data availability and representativeness remain the main challenges for proper determination of Ψmin, a closer look at Ψ distributions and the use of statistically robust meth- ods to estimate Ψmin opens new ground for characterizing plant hydraulic risks. The minimum water potential (Ψmin) describes the largest Introduction xylem tension that a species (or population or individual plant) Extreme values are important in shaping the physiology, ecology experiences in the field. Because in actively transpiring plants and evolution of plants (Gutschick & BassiriRad, 2003). Physio- xylem water potentials decline from roots to leaves along the logical extremes are frequently used to characterize the functional soil–plant–atmosphere continuum, minimum water potentials limits of species and populations, describing the operational lim- are usually determined in leaves or twigs. Low xylem water its of organisms and their tolerance to stressful environmental potentials imperil turgor maintenance and metabolic function factors (e.g. Bozinovic et al., 2011; Kingsolver & Buckley, 2017). (Tyree & Jarvis, 1982; Bartlett et al., 2012) and increase the For this reason, extreme values are frequently used as parameters risk of xylem embolism and the corresponding hydraulic dys- constraining key variables in ecological models, particularly in function (Tyree & Sperry, 1988; Tyree & Zimmermann, the context of predicting responses to ongoing climate change. In 2002). Our current understanding is that hydraulic failure as a plant ecophysiology, examples of such extreme values include result of xylem embolism triggers drought-induced mortality in maximum (Hoshika et al., 2018) and minimum stomatal con- plants (Adams et al., 2017; Choat et al., 2018; Brodribb et al., ductance (Duursma et al., 2019), as well as the minimum tissue 2020). Consequently, the hydraulic safety margin (HSM), usu- water potential characterizing extreme exposure to drought ally quantified as the difference between Ψmin and embolism (Bhaskar & Ackerly, 2006). resistance (e.g. ΨPLC50, or the water potential causing 50% loss 404 New Phytologist (2021) 232: 404–417 Ó 2021 The Authors www.newphytologist.com New Phytologist Ó 2021 New Phytologist Foundation This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made. New Phytologist Methods Research 405 of hydraulic conductivity), has become a standard metric for Clearly, a more robust way of determining Ψmin (and other characterizing plant vulnerability to drought (Meinzer et al., ecophysiological extremes) is needed. Percentiles may be used 2009; Choat et al., 2012; Delzon & Cochard, 2014). This to characterize extreme values from a Ψ distribution, but they variable has the merit of integrating a measure of absolute make use of all observations and are hence highly dependent stress tolerance determined in the laboratory (ΨPLC50) with a on the representativeness of the sample. Extreme value theory measure of (extreme) exposure at the tissue level (Ψmin), thus (EVT) is the branch of statistics devoted to the study of yielding a promising measure of mortality risk (Anderegg extreme deviations from the median, focusing on the tails of et al., 2016; Benito Garzon et al., 2018). Note, however, that the underlying probability distributions as a means to estimate hydraulic failure may be associated with other metrics such as extreme values (Coles, 2001; Reiss & Thomas, 2007; the water potential causing 88% loss of conductivity, which is Embrechts et al., 2013). The two most common approaches frequently used to calculate HSM in angiosperms (Urli et al., involve grouping data into blocks of equivalent sizes and then 2013), or the remaining water transport capacity in the xylem fitting the distribution of block maxima (or minima) using the (Venturas et al., 2016). In addition, leaf shedding may compli- generalized extreme value (GEV) distribution, or using the cate the link between measures of hydraulic risk and mortality, generalized Pareto distribution (GPD) to fit the tail of the dis- particularly in drought-deciduous species (Kolb & Davis, tribution of raw values (‘peak over threshold’ or POT method; 1994). cf. Fig. 1) (see Supporting Information Notes S1 for an The determination of Ψmin is challenging. Some studies extended introduction to EVT methods). have modelled Ψmin from estimated minimum soil water con- Extreme value theory methods are widely applied in finance, tent using pedotransfer functions and information on soil tex- hydrology and climatology, among other fields, but their use in ture (Benito Garzon et al., 2018). This approximation has the ecology is rare and basically limited to the characterization of advantage of providing spatially explicit estimates of Ψmin (and environmental extremes (Katz et al., 2005; Denny et al., 2009) HSM), but it fails to account for the fact that Ψmin is also and the likelihood of species extinction (Burgman et al., 2012). driven by several, covarying plant traits (Bhaskar & Ackerly, We are only aware of two studies that apply EVT to address rele- 2006; Martınez-Vilalta & Garcia-Forner, 2017). A more com- vant questions in plant functional ecology, one to characterize mon approach is to estimate Ψmin as the absolute minimum seed mass distributions (Edwards et al., 2015) and the other to water potential measured in a given species or population. determine maximum stomatal conductance (Murray et al., However, this is also problematic, as sample extremes are 2019). inherently biased and the magnitude of the bias depends on Here, we take advantage of EVT to provide a general frame- the probability distribution of the underlying variable and on work aimed at improving the estimation of Ψmin. We first charac- sample size (Head et al., 2012). In particular, if Ψmin is esti- terize the statistical properties of water potential distributions mated as the absolute minimum Ψ recorded for a given and their variability within species, focusing on distribution tails species, its magnitude will necessarily get more negative (more and taking advantage of a new database of raw water potential extreme) as more Ψ measurements are taken, which seems a values assembled for this study. We then compare different meth- highly undesirable property when different species with ods to estimate Ψmin and the associated uncertainties and biases unequal sample size are compared. Importantly, sampling using two empirical databases and synthetic data characterizing effort is likely to vary directionally with drought exposure typical water potential distributions. Finally, we assess the impli- because species in drier and more variable environments tend cations of our results for the characterization of hydraulic risk to be sampled more intensively than species from wetter envi- and its variability among species, and provide guidelines to esti- ronments, where drought is not perceived as an important mate Ψmin and other ecophysiological extremes in a statistically stress factor. robust manner. An alternative approach to estimating Ψmin involves determin- ing the Ψ associated with mortality. When this is done through Materials and Methods opportunistic sampling under extreme drought conditions in the field, the problem becomes the identification of the physiological Empirical water potential datasets point of no return (cf. Hammond et al., 2019; Sapes et al., 2019), thus avoiding measurements taken in drying tissue once plants The first dataset was collected by Martınez-Vilalta et al. (2014) are effectively dead. This can be achieved by monitoring mea- and hence we identify it as the ‘MV14’ dataset. It consists of a sured branches to ensure they survive the drought (e.g. Venturas global collection of paired predawn and midday leaf water poten- et al., 2016), but this is not always feasible. An alternative tial seasonal time courses for 159 species (3562 individual obser- approach is to conduct drought simulation experiments in con- vations in total). Note that we use the term leaf Ψ for simplicity, trolled environments and apply rewatering trials to identify the but all measurements we refer to in the manuscript were con- point of no return (e.g. Brodribb & Cochard, 2009; Urli et al., ducted in leaves or (leafed) twigs. These data were extracted from 2013). Even these types of study, however, may bring plants published papers and hence correspond to averages across indi- beyond any conditions they could actually experience in the field, viduals for a given site and date, as given typically in figures. This and thus it seems better to refer to the corresponding magnitude dataset was used to study the impact of sample size on the estima- as Ψlethal (cf. Liang et al., 2021). tion of Ψmin for the subset of 17 species with a reasonably large Ó 2021 The Authors New Phytologist (2021) 232: 404–417 New Phytologist Ó 2021 New Phytologist Foundation www.newphytologist.com New 406 Research Methods Phytologist Fig. 1 Time series of measured midday leaf water potentials (top row), histograms showing the probability distributions of measured water potentials (middle row), and tails of the distributions (bottom row) for four representative species included in the ‘raw’ dataset. Each column corresponds to the species indicated on top. The dashed red lines in top and middle rows indicate the 90th percentile of the distribution, which is the threshold value used to define the tail in the bottom panel. The blue line in the bottom panel indicates the generalized Pareto distribution (GPD) fit to the tail of the corresponding distribution. Note the use of log scale in the y-axis of the bottom row panels. Data for Juniperus monosperma and Pinus edulis are pooled across experimental treatments. sample size (N ≥ 50 midday measurements, spanning one to 11 the global xylem functional traits (XFT) dataset by Choat et al. sites depending on the species) (Table S1). (2012) incorporating newer references covering primarily China, A second empirical dataset including raw, individual-level the Amazon and Australia (Sanchez-Martinez et al., 2020). After measurements of midday water potentials was collected specifi- aggregating to the species level this database contained Ψmin data cally for this study, using our own measurements and contributed for 654 species (‘species-level’ dataset). Unlike the measurements or publicly available databases. This dataset contains a total of in the ‘MV14’ and ‘raw’ datasets, which were taken on uncovered 7749 individual measurements of midday leaf water potential on leaves (or twigs), the Ψmin data in the XFT dataset combine mea- 76 species, and we identify it as the ‘raw’ dataset. We selected surements taken on covered and uncovered leaves. However, the data for the 26 species for which the number of measurements effect of the pressure drop in leaves associated with transpiration was ≥ 30 and there were at least 10 combinations of site and sam- is minimal in our case because Ψmin is usually recorded under pling date (e.g. two sites with five or more time points per site or extreme conditions when stomata are closed (cf. Choat et al., one site with 10 or more time points) (Table S2). Note that an 2012; Martin-StPaul et al., 2017). explicit separation of temporal and spatial effects would have been desirable in all analyses, but was not possible owing to data Synthetic datasets limitations. Absolute Ψmin values from the two previous databases We generated two synthetic water potential datasets with realistic (‘MV14’ and ‘raw’) were combined with an updated version of statistical properties to conduct a detailed analysis of the impact New Phytologist (2021) 232: 404–417 Ó 2021 The Authors www.newphytologist.com New Phytologist Ó 2021 New Phytologist Foundation New Phytologist Methods Research 407 of sample size and different analytical approaches on species-level combinations we compared the full Ψ distributions among com- estimates of Ψmin. One of these datasets was generated from a binations of site and year using a k sample Anderson–Darling test Gamma distribution with scale and shape parameters equal to (R package KSAMPLES v.1.2.9), which gives more weight to the those obtained for the fitted Juniperus monosperma distribution, tails of the distributions than do alternative methods (Scholz & and simulated a right-skewed distribution typical of relatively Stephens, 1987). In addition, we also compared the Ψmin values anisohydric species with looser water potential regulation estimated for different combinations of site and year (within (Fig. S1). The second dataset was obtained using a Weibull distri- species) using the following methods: the lowest midday leaf bution with parameters equal to those corresponding to Pinus water potential (highest absolute value) recorded for a given pinaster, and represents a relatively isohydric species with tighter species (cf. Choat et al., 2012); the 99th percentile of the sampled water potential regulation (Fig. S1). values for the species; and the 99th percentiles of the Gamma, Each of the synthetic datasets corresponded to a combination Weibull, and GPD distributions (parameterized as described ear- of data for two to 100 sites, each of which had been measured 10 lier) fitted to the water potential data. times and with 10 plants sampled per time point, and hence pro- vided an ‘ideal’ species-level coverage relative to what is possible Effect of sample size and estimation method on Ψmin determi- in practice in field studies. Individual data points within a sam- nation (species level) To assess the impact of sample size and pling event and seasonal measurements within a year were sam- estimation method on the determination of species-level Ψmin we pled from the same distribution (constant scale and shape used both the ‘MV14’ and ‘raw’ datasets. As earlier, Ψmin was cal- parameters). We also conducted additional simulations in which culated as: the lowest midday leaf water potential, the 99th per- variability was allowed among sites to account for climatic as well centile of the sample, or the 99th percentile of the Gamma, as ecotypic variability. This variability was introduced by sam- Weibull and GPD distributions fitted to the species data. In all pling the corresponding distribution parameters from multivari- cases we obtained Ψmin estimates from the whole dataset for each ate normal distributions with the same means as discussed earlier species and from random subsamples, including 80%, 60%, and a realistic covariance matrix, based on the distribution fits for 40%, 20% and 10% of the data, as well as for N = 10 (in that lat- species with data from multiple populations/years. ter case to compare a common, low sample size across species). Computations were repeated 100 times per species using differ- ent random subsamples. Data analysis To account for all sources of variability in water potential Water potential distributions and their tails (species level) We observations in an idealized situation, we conducted a similar used the ‘raw’ dataset (including completely disaggregated, analysis using the synthetic datasets. In this case the underlying individual-level measurements) to study the tails of Ψ distribu- probability distribution was known and hence the methods to tions (data above the 90th percentile). We compared the fit pro- estimate Ψmin that we compared were: the lowest midday leaf vided by three different probability functions with different tail water potential observed for a given ‘species’; the average Ψmin shapes, corresponding to the log-normal (heavy-tailed), Gamma across populations (cf. Bartlett et al., 2016); the 99.9th percentile and Weibull (heavy-tailed depending on the shape parameter) of the sampled values; the 99.9th percentile of the underlying dis- distributions, with the fit provided by the GPD (using the proba- tribution used to generate the data; the 90th percentile of the bility weighted moment method). We also compared GPD fits GEV distribution fitted to the Ψmin data, blocking by popula- using alternative thresholds to characterize the distribution tail, tion; and the 99.9th percentile of the GPD fitted to the water corresponding to 75th and 95th percentiles. In all analyses, water potential data. Here we used more extreme percentiles than with potential values were converted to positive (absolute) values the empirical datasets (99.9th instead of 99th as a reference) before fitting the different probability distributions. The because sample size was much larger and included all sources of goodness-of-fit of different distributions was compared qualita- variability (across individuals, time and populations). Similarly, tively using the Akaike information criterion (AIC). the 90th percentile was used when fitting the GEV distribution (instead of the 99.9th percentile) because in that case we are fit- Variance components and variability in Ψ and Ψmin within ting (blocked) Ψmin data as opposed to raw Ψ data. Because we species We used variance components analysis to estimate the are blocking by population, each of which corresponds to 100 Ψ contribution of species, site, year and date (day of the year) to values, the 90th percentile of the resulting Ψmin should be compa- total Ψ variability in the ‘raw’ dataset. We conducted additional rable to the 99.9th percentile of the underlying raw data. analyses focusing on the species with largest sample size (J. monosperma and Pinus edulis) to further explore the components Estimating hydraulic risk We used the disaggregated data of of temporal variability. In those models date was nested within the ‘raw’ dataset (full sample size) to calculate the species-level season and season within year. All variance components analyses HSM using different estimates of Ψmin, including the lowest were conducted using the R package VCA (v.1.4.3). Ψmin value, the average Ψmin across populations and years, the To study the variability within species, we subset the data from 99th percentiles of the Weibull, Gamma and GPD distributions, the ‘raw’ dataset by combinations of site and year and selected and the 90th percentile of the GEV distribution fitted to the Ψmin only the combinations with sufficient data (N ≥ 30 individual data, blocking by population and year (the latter only for the only measurements). For each of the six species that had several such two species for which the number of blocks was > 10: Pinus Ó 2021 The Authors New Phytologist (2021) 232: 404–417 New Phytologist Ó 2021 New Phytologist Foundation www.newphytologist.com New 408 Research Methods Phytologist halepensis and Quercus ilex). These estimates were computed for bimodal distributions. This pattern was observed in some species the 17 species for which stem ΨPLC50 data were present in the for which several populations were included (e.g. Quercus coc- ‘species-level’ database. We also used the fitted GPD and GEV cifera), but also in species that had been measured over one single distributions to calculate the probability of Ψmin reaching a value site and year (e.g. Encelia actoni, Ericameria linearifolia and Pur- equal to ΨPLC50. To get a ΨPLC50 value representative of the max- shia tridentata; Fig. S3). As expected, when we focused on the imum resistance to embolism at the species level, and also to tails of the distributions (‘raw’ dataset) the best fit was always pro- avoid methodological issues that may overestimate vulnerability vided by the GPD, which had AIC values that were on average to embolism (Cochard et al., 2013; see also Venturas et al., 13, 14 and 24 units lower than those corresponding to the log- 2019), in these analyses we estimated ΨPLC50 as the minimum normal, Gamma and Weibull distributions, respectively value recorded for a given species. This minimum ΨPLC50 was (Table S3). Fig. 1 shows GPD fits to distribution tails for four tightly related to the species mean ΨPLC50 (r = 0.93, slope = representative species. The average value of the shape parameter 0.85  0.09) and to more extreme measures of hydraulic loss of the GPD distribution was 0.38 and it was < 0 in 19 out of such as ΨPLC88 (r = 0.91, slope = 1.10  0.15), and our results 26 species (Table S3), indicating that the distribution of values remained qualitatively identical if these other metrics were used above the Ψ threshold was bounded in most cases (see Embrechts (data not shown). et al., 2013). However, this shape parameter (unlike the corre- sponding, estimated Ψmin) was relatively unstable and depended Assessing implications for global patterns of hydraulic safety on the choice of threshold (Table S3), probably reflecting sample margins Finally, we assessed whether methodological issues in size limitations. Ψmin estimation could be affecting the global relationship between Ψmin and ΨPLC50 as reported by Choat et al. (2012) and Variance components and within-species variability similar studies. To that end, we first assessed the relationship between ΨPLC50 and sample size (for Ψmin) for both the species The variance components analysis indicated that most of the in our ‘raw’ dataset and globally using the updated XFT database variability in Ψ corresponded to differences among dates (with individual site-specific values retained instead of being (57%) and sites (19%), whereas species accounted for 10% of aggregated by species as in our ‘species-level’ dataset). A positive the total variability (Table S4). A similar pattern emerged relationship between |ΨPLC50| and sample size would indicate within species, where > 50% of the variability was observed that more resistant species are more intensively sampled, and among dates (within seasons and years) for the two species hence a potential, directional bias in the estimation of Ψmin. We tested (Table S4). For five out of six species analysed, the dis- then compared the relationship between Ψmin and ΨPLC50 for the tribution of Ψ differed among years/sites (k sample Anderson– species in our ‘raw’ dataset using two different estimates of Ψmin: Darling test, P < 0.001), the only exception being Phillyrea lat- the absolute minimum (as in Choat et al. (2012)) and the GPD- ifolia, for which we could only compare two consecutive years based estimate described earlier, which is expected to be more at the same site (Fig. 2). For three species (P. halepensis, robust to sample size effects. In both models |Ψmin| and |ΨPLC50| P. sylvestris and Q. ilex) we could also group the observations were log-transformed to normalize the distribution of the residu- by years within a site and compare the distributions by sites. als. In this case ΨPLC50 was estimated as the mean ΨPLC50 for a The resulting Ψ distributions also differed significantly among given species, for consistency with the original analysis by Choat sites (P < 0.001). The fact that Ψ distributions differed within et al. (2012), although other measures of hydraulic vulnerability species also resulted in substantially different Ψmin estimates were also tested (minimum ΨPLC50, ΨPLC88). Note that we could among sites and years, regardless of the estimation method that not perform a similar analysis using the whole ‘species-level’ was used (Table S5). For estimates based on the GPD distri- dataset because the real sample size (number of individuals, sam- bution, which again provided the best fit to the Ψ tail data pling dates and sites used to estimate Ψmin) is unknown in this (Table S5), Ψmin values varied by as much as 5.7 MPa between case. different combinations of site by year in P. halepensis, and up All analyses were conducted in the R environment (R v.3.6.2; to 4.6 MPa when comparing the same population (Orihuela) The R Foundation for Statistical Computing). Distributions between two consecutive years. were fitted using the packages FITDISTRPLUS (v.1.1.1), EVD (GEV, v.2.3.3) and FEXTREMES (GPD, v.3042.82). Comparing Ψmin estimation methods Species-level Ψmin differed substantially across estimation meth- Results ods, both for the ‘raw’ (cf. values corresponding to maximum sample size in Fig. 3) and for the ‘MV14’ datasets (Fig. S4). The statistical distribution of water potential data These values were typically more extreme for either the 99th per- The distribution of leaf water potentials for the species included centile of the Gamma distribution fitted to the data or for the in the ‘MV14’ and ‘raw’ datasets are shown in Figs S2 and S3. As lowest measured Ψmin, with the latter being generally more expected, the shapes of the distributions differed across species extreme for species with a large sample size (N > 100). The Ψmin and were not always well captured by standard distributions (log- estimates based on the 99th percentile of the GPD were usually normal, Gamma, Weibull). In particular, some species showed similar to those obtained from the sample 99th percentile, but the New Phytologist (2021) 232: 404–417 Ó 2021 The Authors www.newphytologist.com New Phytologist Ó 2021 New Phytologist Foundation New Phytologist Methods Research 409 Fig. 2 Density distributions for different combinations of site and year for six representative species from the ‘raw’ dataset. Asterisks at the top left corner of a given panel indicate significant differences among sites/years for the corresponding species (***, P < 0.001). Data for Juniperus monosperma and Pinus edulis are pooled across experimental treatments. former were less sensitive to sample size. The GPD-based esti- Quercus coccifera in Fig. 3; Q. ilex in Fig. S4). The bias associated mates were very robust to the criteria used to set the threshold, with sample size was particularly conspicuous for heavily sampled with resulting Ψmin values differing by 0.05 and 0.03 MPa on species, for which estimates could differ by several MPa (e.g. average between our base calculations (90th percentile) and those between c. 6 and 10 MPa for J. monosperma; Fig. 3). As using the 75th or 95th percentiles, respectively (Table S3). In expected, Ψmin estimates based on the 99th percentile of the sam- addition, GPD-based estimates were also reasonably robust to pled distribution were basically identical to the absolute mini- removing different proportions of the data with less extreme Ψ mums for sample sizes below 100, and shared the strong values (i.e. using only the data measured under relatively extreme dependency on sample size. All the other Ψmin estimates, based conditions) (Table S6). on percentiles of the fitted Gamma, Weibull or GPD distribu- tions, were much less sensitive to sample size. However, they were all affected by sample size at low values (N < 20), at least for some Effect of sample size on Ψmin estimates species (Figs 3, S4). Sample size had a strong effect on Ψmin when it was estimated as Sample size also had a strong effect on Ψmin estimated as the lowest midday leaf water potential recorded for a given the lowest midday leaf water potential recorded in the syn- species, regardless of whether data came from the ‘raw’ or the thetic datasets, regardless of the distribution used to generate ‘MV14’ dataset (Figs 3 and S4, respectively). In most species the the data (Fig. 4). At high sample sizes, these estimates were effect tended to saturate at high sample sizes, but in some cases a always substantially more extreme than those obtained by the plateau was not observed even for sample sizes > 100 (e.g. other methods, including the 99.9th percentile from the known Ó 2021 The Authors New Phytologist (2021) 232: 404–417 New Phytologist Ó 2021 New Phytologist Foundation www.newphytologist.com New 410 Research Methods Phytologist Fig. 3 Effect of sample size on minimum water potential estimates (absolute values in MPa) from the ‘raw’ dataset. For each species, different estimation methods are compared and shown in different colours (black lines: the absolute lowest value from the sample; grey lines: the 99th percentile from the sample, which is indistinguishable from absolute lowest values for sample sizes < 100; red lines: the 99th percentiles of the Gamma distribution fitted to the data; blue lines: the 99th percentiles of the Weibull distribution fitted to the data; green lines: the 99th percentile of the generalized Pareto distribution (GPD) fitted to the data). Error bars show SDs of 100 random repetitions for each species and sample size. distribution. At the other extreme, when ‘species’ Ψmin was Estimating hydraulic risk estimated as the average Ψmin across populations, the resulting values were very stable but always substantially less extreme The calculated HSMs for the 17 species in the ‘raw’ dataset for than the 99.9th percentile from the known distribution which stem ΨPLC50 was available were highly variable depending (Fig. 4). Estimates obtained from the sample 99.9th percentile on the method used to estimate Ψmin, with discrepancies of more and those using the GEV and GPD distributions all converged than 1 MPa in many species and up to > 4 MPa in P. halepensis to values indistinguishable from the 99.9th percentile of the (Fig. 5). For all the species in which there was more than one original distribution (Fig. 4). These three methods were quite combination of year and site, the HSM estimates using the abso- robust to low sample sizes down to N = 10 populations, lute lowest Ψmin were substantially lower than those using the although the estimate from the sample 99.9th percentile typi- mean Ψmin across sites/years. The HSM estimates obtained using cally showed the highest standard deviations. For even lower the Weibull and particularly the Gamma distributions tended to sample sizes, GPD-based estimates appeared the most consis- be substantially lower than those using the GPD distribution, tent, although uncertainty increased substantially in all cases, whereas those based on the GPD and GEV distributions were particularly for GEV-based estimates. These results remained similar for the two species for which the latter could be computed qualitatively similar when variability in the parameters of the (Fig. 5). underlying distributions was introduced among populations The probability of individual water potential measurements (Fig. S5). reaching ΨPLC50 calculated from the GPD (90th percentile New Phytologist (2021) 232: 404–417 Ó 2021 The Authors www.newphytologist.com New Phytologist Ó 2021 New Phytologist Foundation New Phytologist Methods Research 411 Fig. 4 Effect of sample size (number of sites, from two to 100) on minimum water potential estimates (absolute values in MPa) from the synthetic datasets. Each site comprised data for 10 measurement dates with 10 plants sampled per time point. Results are shown for two different distributions: a Gamma (left panel) and a Weibull distribution (right panel) (see text for details on parameterization and Supporting Information Fig. S1 for the corresponding histograms). Different estimation methods are compared and shown in different colours (red: the absolute lowest value; brown: the mean of the minimum values observed at different sites; blue: the 99.9th percentile from the whole sample; black: the 99.9th percentile of the known distribution; green: the 99.9th percentile of the generalized Pareto distribution (GPD) fitted to the data; orange: the 90th percentile of the generalized extreme value (GEV) distribution fitted to the data blocked by site). Error bars show SDs of 100 random repetitions for each species and sample size. Note that some error bars for the GEV estimates under low sample sizes were very large and have been omitted. threshold) ranged between 0% for four species to > 10% for if water potential gradients in the xylem are expected to be low Ambrosia salsola, Rhus ovata, Viburnum tinus and Vitis vinifera under extreme drought (Choat et al., 2012). (Table S7). The latter values imply that ΨPLC50 was higher (less negative) than the threshold used to fit the GPD, and hence the Implications for global patterns of hydraulic safety margins probability of reaching ΨPLC50 could not be estimated with preci- sion. The consistency of these probability estimates when com- We found a significant, negative correlation between species-level paring different GPD thresholds (75th and 95th percentiles) was ΨPLC50 and the number of entries in the ‘species-level’ dataset (r quite variable across species (4.9–150% deviation depending on = 0.22, P < 0.001, N = 226), suggesting that embolism- the species and threshold comparison; Table S7), and was clearly resistant species tended to be more intensively sampled. A similar lower than for the Ψmin estimates themselves (Table S3). For the relationship was found for the species in our ‘raw’ dataset, which two species for which GEV distributions could be fitted, the was statistically significant when minimum ΨPLC50 was consid- probability of Ψmin reaching ΨPLC50 (or lower) differed substan- ered (r = 0.55, P = 0.02, N = 17) but not when using mean tially (3% for Q. ilex, 19% for P. halepensis). This difference was ΨPLC50 (r = 0.35, P = 0.17, N = 17). Consistent with the origi- qualitatively consistent with that obtained using GPD-based esti- nal analysis in Choat et al. (2012) we found a strong positive rela- mates (Table S7) and reflected differences in the corresponding tionship between Ψmin and mean ΨPLC50, regardless of whether HSM (Fig. 5). Ψmin was estimated from absolute minimum values or GPD- We note that our focus here is on the effect on HSM of differ- based estimates (r = 0.65, P < 0.01, N = 17 in both cases) (Fig. 6). ent approaches to estimate Ψmin, not on the HSM values them- The slope of the relationship was significantly < 1 in the two selves, which should be interpreted with caution. In particular, cases, but declined from 0.44  0.13 (intercept = 0.94  0.19) HSM estimates for some species may be problematic owing to when Ψmin was estimated from absolute minimum values to uncertainty in ΨPLC50 values. For instance, the highly negative 0.38  0.12 (intercept = 0.92  0.16) when GPD-based esti- HSM in Viburnum tinus, and particularly Vitis vinifera, probably mates were used (Fig. 6). These results indicate that more resis- arise from the overestimation of embolism vulnerability for these tant species had wider safety margins, and that this trend was species (Charrier et al., 2018), for which only one ΨPLC50 value exacerbated when more robust Ψmin estimates that accounted for was available in the XFT dataset. The fact that we are comparing sample size effects were used. We fitted a single relationship for leaf/twig Ψmin with stem Ψ50 may introduce additional bias, even angiosperms and gymnosperms because, in our dataset, HSM Ó 2021 The Authors New Phytologist (2021) 232: 404–417 New Phytologist Ó 2021 New Phytologist Foundation www.newphytologist.com New 412 Research Methods Phytologist Fig. 5 Comparison of hydraulic safety margins (HSM, MPa) obtained using different estimates of minimum water potential (Ψmin) on the species from the ‘raw’ dataset. Different estimation methods include the absolute lowest value from the sample (‘Abs’), the mean across all the sites/years sampled (‘Mean’), the 99th percentile from the fitted Gamma distribution (‘Gamma p99’), the 99th percentile from the fitted Weibull distribution (‘Weibull p99’), the 99th percentile of the fitted generalized Pareto distribution (‘GPD p99’), and the 90th percentile of the generalized extreme value distribution fitted to site/year blocks, for the two species for which the number of combinations of site and year was > 10 (‘GEV p90’). HSM is quantified as the difference between Ψmin and stem ΨPLC50, and hence a lower (e.g. negative) value indicates higher hydraulic risk. Please note that HSM values reported in this figure are primarily intended for comparison of alternative Ψmin estimation methods and should be interpreted with caution (see text). were similar between these two group of species (P = 0.23, magnitude is generally unknown. Moreover, this bias is not ran- N = 11 angiosperm and N = 6 gymnosperm species). domly distributed across species. Because species living in drier environments tend to be more intensively sampled, current Ψmin estimates for these species would be necessarily more extreme Discussion (even if the underlying probability distribution was identical) than those for species living under milder conditions. Current estimates of Ψmin are biased The fact that Ψmin estimates are biased has significant implica- Minimum water potential (Ψmin) is a key variable to characterize tions for the determination of HSMs. Because HSM is usually dehydration tolerance and hydraulic risk in plants (Choat et al., calculated as the difference between Ψmin and a measure of 2012, 2018; Delzon & Cochard, 2014). Our results show that hydraulic resistance (e.g. ΨPLC50), estimated HSM would have current Ψmin values, estimated as the absolute minimum leaf the same absolute bias as the corresponding Ψmin values (cf. water potential reported for a given species, are biased. This is Fig. 5). We also show that this bias can have important conse- because absolute extremes are a function of sample size and the quences for our understanding of the distribution of drought underlying statistical distribution, which differ markedly among resistance and hydraulic risks at the global scale. Choat et al. species (and among sites and years within species; Fig. 2). This (2012) showed that, for angiosperms, HSM were relatively bias can be large, of the order of several MPa (Figs 3, 5), and its invariant among biomes and only increased slightly for more New Phytologist (2021) 232: 404–417 Ó 2021 The Authors www.newphytologist.com New Phytologist Ó 2021 New Phytologist Foundation New Phytologist Methods Research 413 (Breshears et al., 2008; McDowell et al., 2008; Meinzer et al., 2014). The Ψ distribution, however, is also determined by envi- ronmental conditions (Hochberg et al., 2018), which probably explains part of the within-species variability we observed (Fig. 2). This effect could be accounted for by modelling leaf from soil Ψ (e.g. Martınez-Vilalta et al., 2014) or by directly focusing on the difference between the two. This type of analysis should be facilitated by increased availability of environmental drivers and detailed soil and plant water content measures from remote sensing (cf. Wu et al., 2021). Guidelines for the estimation of Ψmin and other ecophysiological extremes Fig. 6 Relationship between water potential at 50% loss of hydraulic Any proper estimation of extreme values from field measure- conductivity (ΨPLC50) and minimum water potential (Ψmin) across the ments has to deal with the fact that these values are rare by species in the ‘raw’ dataset. Red dots correspond to the absolute lowest definition, and thus their estimation is less straightforward than Ψmin recorded for a given species, whereas green dots show generalized Pareto distribution (GPD)-based estimates of Ψmin, which are robust to for typical measures of centrality (e.g. the mean) or dispersion sample size effects. Observations for a given species are connected by a (e.g. the standard deviation). Although the absolute minimum black arrow. Different symbols are used for angiosperm and gymnosperm of any given sample is easy to calculate, this magnitude does species. Red and green lines show the corresponding linear fits. The 1 : 1 not provide a sound basis to estimate Ψmin, as it continuously line is indicated by the dotted black line. Note the log scales in both axes. declines (becomes more extreme) as sample size increases, and hence cannot be used as a reference to compare species (or embolism-resistant species. Our analyses show that this interpre- populations) in practice. An operational definition of Ψmin is tation may need to be reassessed, as HSM estimates became wider therefore needed, and such definition should be probabilistic when more robust, GPD-based estimates of Ψmin were used, rela- to make it independent of sampling effort. Although directly tive to using absolute lowest Ψ values, and this happened particu- estimating (and comparing) a given sample percentile could be larly in more embolism-resistant species as a result of their higher a reasonable solution in many cases, these estimates are more sampling effort (Fig. 6). Although it is impossible to quantify this sensitive to low sample sizes (Fig. 4), and would be particularly effect precisely at present, because the sample size used to esti- problematic when the underlying probability distribution is mate Ψmin in global hydraulic databases is unknown, these results heavy-tailed. We thus advocate here for the use of EVT (Coles, suggest that safety margins may actually increase substantially 2001; Reiss & Thomas, 2007) as a good basis for improved with embolism resistance at the global scale. This is consistent estimation. with the higher difference between the water potential at stomatal Extreme value theory focuses on the tail of the probability dis- closure and ΨPLC50 (another measure of the safety margin) in tribution, as the tail conveys the most relevant information to embolism-resistant species (Martınez-Vilalta et al., 2014; Men- estimate the magnitude and likelihood of extreme events. In gen- cuccini et al., 2015; Martin-StPaul et al., 2017). We note, how- eral, the most efficient approach uses the ‘peak over threshold’ ever, that these patterns are likely to be scale-dependent, as (POT) method (Coles, 2001). This approach involves identifying increased hydraulic risk in embolism-resistant species has been a threshold in the raw distribution of Ψ that characterizes extreme reported in local communities (e.g. Venturas et al., 2016). situations, fitting a GPD to values exceeding that threshold (the The fact that the distribution of leaf water potentials differs tail of the distribution), and then estimating the species’ Ψmin as greatly among species has been largely overlooked in the litera- a given (conventional) percentile of the fitted GPD. Because this ture. Besides having substantial implications for the estimation of method focuses on the most extreme values of the overall distri- Ψmin, this variation reflects important biological phenomena. In bution, it has the advantage of being relatively insensitive to vari- particular, the distribution of leaf water potentials is the end ability in Ψ among populations or years (cf. Fig. 2), which is result of all the control mechanisms by which plants regulate particularly useful when the objective is to obtain species-level their water status as a response to environmental drivers, and estimates of Ψmin. We suggest that a 99.9th percentile would be a hence provides a synthetic, statistical description of water poten- reasonable choice, as it corresponds to the minimum value tial regulation (Martınez-Vilalta & Garcia-Forner, 2017). Under recorded, on average, every 1000 individual water potential mea- the same environmental conditions, relatively anisohydric species surements. Based on our analyses, Ψmin estimates obtained this would have Ψ distributions that tend to be skewed to the right way are consistent when the number of individual Ψ measure- and show fatter tails, whereas more isohydric species would be ments exceeds 20 (Fig. 3), and are robust to different threshold characterized by more left-skewed distributions. Examples of this choices (see Coles (2001) for additional criteria to set the GPD would be J. monosperma and P. edulis (or even more clearly threshold). P. pinaster; Fig. S3), which are paradigmatic examples of anisohy- Using the GPD approach, however, requires fitting raw Ψ dric and isohydric regulation of water potentials, respectively data, which are not always available. In those cases, and Ó 2021 The Authors New Phytologist (2021) 232: 404–417 New Phytologist Ó 2021 New Phytologist Foundation www.newphytologist.com New 414 Research Methods Phytologist particularly when midday Ψ measurements are available for sev- < Ψmin except for J. monosperma. The largest difference between eral sites and years with a reasonably balanced experimental Ψlethal and Ψmin was observed for P. edulis (2.12 MPa), for which design, an alternative way to estimate Ψmin involves grouping the the absolute lowest Ψ after c. 1500 measurements (including an data into blocks corresponding to combinations of site and year. extreme drought treatment) was 4.19 MPa, still a long way Then, the absolute minimum Ψ is obtained for each block (anal- from the estimated Ψlethal ( 5.53 MPa on average). This discrep- ogous to data binning procedures), the GEV distribution is fitted ancy illustrates the fact that at least some species may effectively to these block extremes (cf. Makkonen & Tikanm€aki (2019) for never reach, in the field, the Ψlethal values determined under con- a comparison of fitting methods), and species’ Ψmin is estimated trolled experimental conditions. Tagging and monitoring as a given (conventional) percentile of the fitted distribution. In branches sampled under extreme drought to ensure they were this case, the 90th percentile, corresponding to the minimum alive and survived is a promising avenue to reconcile Ψlethal and water potential that would be observed, on average, once every Ψmin estimates. 10 yr (or 10 combinations of site and year), may be adequate in most cases. That implies, of course, that robust estimates of Ψmin Towards a more robust determination of hydraulic risk can only be obtained if data are available for at least 10 combina- tions of site and year (cf. Fig. 4), which unfortunately is quite rare Characterizing the statistical distribution of water potentials and at present (e.g. only two out of the 76 species recorded in our their minimum values opens new ground for the determination ‘raw’ dataset fulfilled this requirement). The GEV-based esti- of plant hydraulic risk, as it allows the probability of reaching a mates may be particularly useful to capture biologically relevant certain, critical Ψ to be estimated. If enough data are available, differences among sites or years. It is noteworthy that, in our case, deriving return intervals analogous to those used in hydrology the resulting Ψmin estimates were similar for GEV- and GPD- and other applied sciences (i.e. the average time between events based methods for the two species in which both could be applied of a given magnitude) becomes straightforward (e.g. Reiss & (Fig. 5). Thomas, 2007). Because hydraulic vulnerability thresholds are Sampling issues and representativeness remain the main diffi- relatively well understood, at least for the xylem, estimating the culty in obtaining robust estimates of Ψmin. Water potential mea- probability of a given amount of hydraulic damage (e.g. 50% or surements tend to be heavily clustered around specific sites or 88% loss of conductivity), or the corresponding return interval, years, and hence may not capture the relevant conditions (severe also becomes possible. Conceptually, this is clearly superior to drought) and spatiotemporal variability to accurately determine HSMs estimated as the difference between Ψmin and embolism extreme values. Our variance components analyses show that resistance, as currently done. Unfortunately, our results show that most of the Ψ variability occurs over time, and particularly at the the quality of species-level, water potential data is rarely good daily scale (Table S4), and hence temporal replication is critical enough at present to obtain robust probabilistic estimates of to characterize Ψ distributions. At the same time, it is clear that hydraulic risk. Therefore, even without considering the issues ecotypical variation may also be high (e.g. Atzmon et al., 2004; associated with the determination of hydraulic vulnerability, see also our Fig. 2) and hence a proper characterization of Ψmin at which are beyond the scope of this paper, a proper characteriza- the species level requires some degree of spatial replication. tion of hydraulic risk remains a challenge (Martinez-Vilalta et al., Because water potential measurements are largely conducted to 2019). study vegetation responses to drought, sampling schemes tend to In conclusion, we demonstrate that current Ψmin estimates prioritize dry locations and dry periods, effectively biasing water are biased and propose two EVT-based methods to obtain sta- potential distributions towards low values. Although, by focusing tistically robust estimates of Ψmin. We validate these methods, on the tail of the Ψ distribution, EVT-based methods are rela- which may also be applied to other ecophysiological extremes, tively insensitive to the impact of undersampling other parts of using different empirical and synthetic datasets at different the distribution, they are not immune to bias, and extensions measurement levels. Our results show the potential of EVT- may be required to account for relevant sources of variation and based methods to improve the determination of hydraulic risks to obtain accurate estimates. Psychrometers offer a powerful in plants, and also point towards novel ways of characterizing alternative to obtain quasi-continuous Ψ series and increase data water potential regulation. At the same time, our analyses representativeness, at least in the temporal dimension (e.g. Guo highlight the large spatiotemporal variability in Ψ within et al., 2020), and have the additional advantage of measuring species, particularly at daily scales, and the limitations of cur- stem (as opposed to leaf) water potential, which is more directly rently available data. We urge improved sampling designs to comparable with most assessments of hydraulic vulnerability. better characterize water potential variability within species and However, psychrometer measurements are still rare and cannot also better reporting of these data, if possible including be obtained for all taxa. In summary, even though sampling individual-level measurements or at least a clear indication of issues are minimized by the EVT approach we propose, they are the level of aggregation. Making raw data available together certainly not solved. with relevant metadata (e.g. site coordinates, time stamps, local Interestingly, for the seven species in the ‘raw’ dataset for environmental data, occurrence of severe drought effects) which we could obtain Ψlethal values from the literature would pave the way for building global, harmonized water (Table S8), those compared well with our Ψmin estimates potential databases and for a rigorous treatment of inclusion (r = 0.76, average difference < 0.9 MPa). For these species, Ψlethal criteria for the determination of Ψmin and associated measures New Phytologist (2021) 232: 404–417 Ó 2021 The Authors www.newphytologist.com New Phytologist Ó 2021 New Phytologist Foundation New Phytologist Methods Research 415 of hydraulic risk. Altogether, this effort should result in better Bartlett MK, Klein T, Jansen S, Choat B, Sack L. 2016. The correlations and understanding and prediction of plant responses to drought sequence of plant stomatal, hydraulic, and wilting responses to drought. Proceedings of the National Academy of Sciences, USA 113: 13098–13103. under ongoing climate change (Choat et al., 2018; Brodribb Bartlett MK, Scoffoni C, Sack L. 2012. The determinants of leaf turgor loss et al., 2020; Ruiz-Benito et al., 2020). point and prediction of drought tolerance of species and biomes: a global meta- analysis. Ecology Letters 15: 393–405. Benito Garzon MB, Mu~ noz NG, Wigneron J-P, Moisy C, Ferna ndez-Manjarres Acknowledgements J, Delzon S, Jordan G. 2018. The legacy of water deficit on populations having experienced negative hydraulic safety margin. Global Ecology and Biogeography This research was supported by the Spanish Ministry of Science 27: 346–356. and Innovation (MICINN) via competitive grants CGL2017- Bhaskar R, Ackerly DD. 2006. Ecological relevance of minimum seasonal water 89149-C2-1-R and CGL2017-89149-C2-2-R. AV thanks potentials. Physiologia Plantarum 127: 353–359. MICINN for the project Inertia (PID201-111332RB-C22). JM- Bozinovic F, Calosi P, Spicer JI. 2011. Physiological correlates of geographic V benefited from an ICREA Academia award. We want to thank range in animals. Annual Review of Ecology, Evolution, and Systematics 42: 155– 179. Sanna Sevanto, Nate McDowell and co-authors for making the Breshears DD, Myers OB, Meyer CW, Barnes FJ, Zou CB, Allen CD, SUMO water potential data publicly available (https://doi.org/ McDowell NG, Pockman WT. 2008. Tree die-off in response to global 10.15485/1439886). change-type drought: mortality insights from a decade of plant water potential measurements. Frontiers in Ecology and the Environment 7: 185–189. Brodribb TJ, Cochard H. 2009. Hydraulic failure defines the recovery and point Author contributions of death in water-stressed conifers. Plant Physiology 149: 575–584. Brodribb TJ, Powers J, Cochard H, Choat B. 2020. Hanging by a thread? JM-V, MM, LS, RP and MdC conceived the ideas and designed Forests and drought. Science 368: 261–266. methodology; RP, LS, IA, SD, AV and JM-V provided data; JM- Burgman M, Franklin J, Hayes KR, Hosack GR, Peters GW, Sisson SA. 2012. V and LB analysed the data; JM-V led the writing of the Modeling extreme risks in ecology. Risk Analysis 32: 1956–1966. manuscript. All authors contributed critically to the drafts and Charrier G, Delzon S, Domec J-C, Zhang Li, Delmas CEL, Merlin I, Corso D, King A, Ojeda H, Ollat N et al. 2018. Drought will not leave your glass gave final approval for publication. empty: Low risk of hydraulic failure revealed by long-term drought observations in world’s top wine regions. Science Advances 4: eaao6969. Choat B, Brodribb TJ, Brodersen CR, Duursma RA, Lopez R, Medlyn BE. ORCID 2018. Triggers of tree mortality under drought. Nature 558: 531–539. Ismael Aranda https://orcid.org/0000-0001-9086-7940 Choat B, Jansen S, Brodribb TJ, Cochard H, Delzon S, Bhaskar R, Bucci SJ, Feild TS, Gleason SM, Hacke UG et al. 2012. Global convergence in the Llorencß Badiella https://orcid.org/0000-0002-9653-7421 vulnerability of forests to drought. Nature 491: 752–755. Miquel Caceres de https://orcid.org/0000-0001-7132-2080 Cochard H, Badel E, Herbette S, Delzon S, Choat B, Jansen S. 2013. Methods Sylvain Delzon https://orcid.org/0000-0003-3442-1711 for measuring plant vulnerability to cavitation: a critical review. Journal of Jordi Martınez-Vilalta https://orcid.org/0000-0002-2332- Experimental Botany 64: 4779–4791. 7298 Coles S. 2001. An introduction to statistical modeling of extreme values. London, UK: Springer-Verlag. Maurizio Mencuccini https://orcid.org/0000-0003-0840- Delzon S, Cochard H. 2014. Recent advances in tree hydraulics highlight the 1477 ecological significance of the hydraulic safety margin. New Phytologist 203: Rafael Poyatos https://orcid.org/0000-0003-0521-2523 355–358. Louis S. Santiago https://orcid.org/0000-0001-5994-6122 Denny MW, Hunt LJH, Miller LP, Harley CDG. 2009. On the prediction of Alberto Vilagrosa https://orcid.org/0000-0002-1432-1214 extreme ecological events. Ecological Monographs 79: 397–421. Duursma RA, Blackman CJ, Lopez R, Martin-StPaul NK, Cochard H, Medlyn BE. 2019. On the minimum leaf conductance: its role in models of plant water Data availability use, and ecological and environmental controls. New Phytologist 221: 693–705. Edwards W, Moles AT, Chong C. 2015. Generalised extreme value distributions All the data sources that support the findings of this study are provide a natural hypothesis for the shape of seed mass distributions. PLoS available in the Supporting Information of this article. The full ONE 10: e0121724. Embrechts P, Kl€ uppelberg C, Mikosch T. 2013. Modelling extremal events for databases are available from the corresponding author upon rea- insurance and finance. Berlin, Germany: Springer-Verlag. sonable request. Guo JS, Hultine KR, Koch GW, Kropp H, Ogle K. 2020. Temporal shifts in iso/anisohydry revealed from daily observations of plant water potential in a dominant desert shrub. New Phytologist 225: 713–726. References Gutschick VP, BassiriRad H. 2003. Extreme events as shaping physiology, Adams HD, Zeppel MJB, Anderegg WRL, Hartmann H, Landh€a usser SM, ecology, and evolution of plants: toward a unified definition and evaluation of Tissue DT, Huxman TE, Hudson PJ, Franz TE, Allen CD et al. 2017. A their consequences. New Phytologist 160: 21–42. multi-species synthesis of physiological mechanisms in drought-induced tree Hammond WM, Yu K, Wilson LA, Will RE, Anderegg WRL, Adams HD. mortality. Nature Ecology & Evolution 1: 1285. 2019. Dead or dying? Quantifying the point of no return from hydraulic Anderegg WRL, Klein T, Bartlett M, Sack L, Pellegrini AFA, Choat B, Jansen failure in drought-induced tree mortality. New Phytologist 223: 1834–1843. S. 2016. Meta-analysis reveals that hydraulic traits explain cross-species Head AW, Hardin JS, Adolph SC. 2012. Methods for estimating peak patterns of drought-induced tree mortality across the globe. Proceedings of the physiological performance and correlating performance measures. National Academy of Sciences, USA 113: 5024–5029. Environmental and Ecological Statistics 19: 127–137. Atzmon N, Moshe Y, Schiller G. 2004. Ecophysiological response to severe Hochberg U, Rockwell FE, Holbrook NM, Cochard H. 2018. Iso/Anisohydry: a drought in Pinus halepensis Mill. trees of two provenances. Plant Ecology 171: plant-environment interaction rather than a simple hydraulic trait. Trends in 15–22. Plant Science 23: 112–120. Ó 2021 The Authors New Phytologist (2021) 232: 404–417 New Phytologist Ó 2021 New Phytologist Foundation www.newphytologist.com New 416 Research Methods Phytologist Hoshika Y, Osada Y, de Marco A, Pe~ nuelas J, Paoletti E. 2018. Global diurnal Tyree MT, Zimmermann MH. 2002. Xylem structure and the ascent of Sap. and nocturnal parameters of stomatal conductance in woody plants and major Berlin, Germany: Springer-Verlag. crops. Global Ecology and Biogeography 27: 257–275. Urli M, Porte AJ, Cochard H, Guengant Y, Burlett R, Delzon S. 2013. Xylem Katz RW, Brush GS, Parlange MB. 2005. Statistics of extremes: modeling embolism threshold for catastrophic hydraulic failure in angiosperm trees. Tree ecological disturbances. Ecology 86: 1124–1134. Physiology 33: 672–683. Kingsolver JG, Buckley LB. 2017. Quantifying thermal extremes and biological Venturas MD, MacKinnon ED, Dario HL, Jacobsen AL, Pratt RB, Davis SD. variation to predict evolutionary responses to changing climate. Philosophical 2016. Chaparral shrub hydraulic traits, size, and life history types relate to Transactions of the Royal Society B: Biological Sciences 372: 20160147. species mortality during California’s historic drought of 2014. PLoS ONE 11: Kolb KJ, Davis SD. 1994. Drought tolerance and xylem embolism in co- e0159145. occurring species of coastal sage and chaparral. Ecology 75: 648–659. Venturas MD, Pratt RB, Jacobsen AL, Castro V, Fickle JC, Hacke UG. 2019. Liang X, Ye Q, Liu H, Brodribb TJ. 2021. Wood density predicts mortality Direct comparison of four methods to construct xylem vulnerability curves: threshold for diverse trees. New Phytologist 229: 3053–3057. Differences among techniques are linked to vessel network characteristics. Makkonen L, Tikanm€aki M. 2019. An improved method of extreme value Plant, Cell & Environment 42: 2422–2436. analysis. Journal of Hydrology X 2: 100012. Wu G, Guan K, Li Y, Novick KA, Feng X, McDowell NG, Konings AG, Martinez-Vilalta J, Anderegg WRL, Sapes G, Sala A. 2019. Greater focus on Thompson SE, Kimball JS, De Kauwe MG et al. 2021. Interannual variability water pools may improve our ability to understand and anticipate drought- of ecosystem iso/anisohydry is regulated by environmental dryness. New induced mortality in plants. New Phytologist 223: 22–32. Phytologist 229: 2562–2575. Martınez-Vilalta J, Garcia-Forner N. 2017. Water potential regulation, stomatal behaviour and hydraulic transport under drought: deconstructing the iso/ anisohydric concept. Plant, Cell & Environment 40: 962–976. Supporting Information Martınez-Vilalta J, Poyatos R, Aguade D, Retana J, Mencuccini M. 2014. A new look at water transport regulation in plants. New Phytologist 204: 105– Additional Supporting Information may be found online in the 115. Supporting Information section at the end of the article. Martin-StPaul N, Delzon S, Cochard H. 2017. Plant resistance to drought depends on timely stomatal closure. Ecology Letters 20: 1437–1447. Fig. S1 Histograms showing the water potential distribution for McDowell N, Pockman WT, Allen CD, Breshears DD, Cobb N, Kolb T, Plaut the two ‘species’ included in the synthetic datasets. J, Sperry J, West A, Williams DG et al. 2008. Mechanisms of plant survival and mortality during drought: why do some plants survive while others succumb to drought? New Phytologist 178: 719–739. Fig. S2 Histograms showing the distribution of water potentials Meinzer FC, Johnson DM, Lachenbruch B, McCulloh KA, Woodruff DR. (absolute value in MPa) from the 17 species in the ‘MV14’ 2009. Xylem hydraulic safety margins in woody plants: coordination of dataset. stomatal control of xylem tension with hydraulic capacitance. Functional Ecology 23: 922–930. Meinzer FC, Woodruff DR, Marias DE, Mcculloh KA, Sevanto S. 2014. Fig. S3 Histograms showing the distribution of water potentials Dynamics of leaf water relations components in co-occurring iso- and (absolute value in MPa) from the 26 species in the ‘raw’ dataset. anisohydric conifer species. Plant, Cell & Environment 37: 2577–2586. Mencuccini M, Minunno F, Salmon Y, Martınez-Vilalta J, H€oltt€a T. 2015. Fig. S4 Effect of sample size on minimum water potential esti- Coordination of physiological traits involved in drought-induced mortality of mates from the ‘MV14’ dataset. woody plants. New Phytologist 208: 396–409. Murray M, Soh WK, Yiotis C, Batke S, Parnell AC, Spicer RA, Lawson T, Caballero R, Wright IJ, Purcell C et al. 2019. Convergence in maximum Fig. S5 Effect of sample size on minimum water potential esti- stomatal conductance of C3 woody angiosperms in natural ecosystems across mates from the synthetic datasets, allowing variability in the bioclimatic zones. Frontiers in Plant Science 10: 558. underlying distribution between sites. Reiss R-D, Thomas M. 2007. Statistical analysis of extreme values, with applications to insurance, finance, hydrology and other fields. Basel, Switzerland: Birkh€auser. Notes S1 A very short introduction to the statistics of extreme Ruiz-Benito P, Vacchiano G, Lines ER, Reyer CPO, Ratcliffe S, Morin X, values. Hartig F, M€a kel€a A, Yousefpour R, Chaves JE et al. 2020. Available and missing data to model impact of climate change on European forests. Ecological Table S1 Species selected from the ‘MV14’ dataset. Modelling 416: 108870. Sanchez-Martinez P, Martınez-Vilalta J, Dexter KG, Segovia RA, Mencuccini M. 2020. Adaptation and coordinated evolution of plant hydraulic traits. Table S2 Species selected from the ‘raw’ dataset. Ecology Letters 23: 1599–1610. Sapes G, Roskilly B, Dobrowski S, Maneta M, Anderegg WRL, Martinez- Table S3 Comparison of different distribution fits to the tails of Vilalta J, Sala A. 2019. Plant water content integrates hydraulics and carbon the water potential distribution data for all the species in the depletion to predict drought-induced seedling mortality. Tree Physiology 39: ‘raw’ dataset. 1300–1312. Scholz FW, Stephens MA. 1987. K-sample anderson-darling tests. Journal of the American Statistical Association 82: 918–924. Table S4 Variance components analysis of Ψ for all the species in Tyree MT, Jarvis PG. 1982. Water in tissues and cells. In: Lange OL, Nobel PS, the ‘raw’ dataset and specifically for Juniperus monosperma and Osmond CB, Ziegler H, eds. Encyclopedia of plant physiology. Physiological plant Pinus edulis. ecology II: Water relations and carbon assimilation. Berlin/Heidelberg, Germany: Springer, 35–77. Tyree MT, Sperry JS. 1988. Do woody plants operate near the point of Table S5 Comparison of different distribution fits for different catastrophic xylem dysfunction caused by dynamic water stress? Answers from a combinations of site and year for six representative species from model. Plant Physiology 88: 574–580. the ‘raw’ dataset. New Phytologist (2021) 232: 404–417 Ó 2021 The Authors www.newphytologist.com New Phytologist Ó 2021 New Phytologist Foundation New Phytologist Methods Research 417 Table S6 Comparison of Ψmin values estimated using the GPD Table S8 Comparison of Ψmin values estimated in this study with approach on different data subsets for Juniperus monosperma and Ψlethal values obtained from literature sources. Pinus edulis. Please note: Wiley Blackwell are not responsible for the Table S7 Mean water potential causing 50% loss of hydraulic content or functionality of any Supporting Information conductivity (ΨPLC50), minimum ΨPLC50, mean water potential supplied by the authors. Any queries (other than missing causing 88% loss of hydraulic conductivity (ΨPLC88) and proba- material) should be directed to the New Phytologist Central bility of Ψmin reaching ΨPLC50 for the species in the ‘raw’ dataset. Office. Ó 2021 The Authors New Phytologist (2021) 232: 404–417 New Phytologist Ó 2021 New Phytologist Foundation www.newphytologist.com