This is the accepted version of the journal article: Sanchez Martinez, Pablo; Martínez Vilalta, Jordi; Dexter, Kyle G.; [et al.]. «Adaptation and coordinated evolution of plant hydraulic traits». Ecology Letters, Vol. 23, Issue 11 (November 2020), p. 1599-1610. DOI 10.1111/ele.13584 This version is available at https://ddd.uab.cat/record/257069 under the terms of the license TITLE: Adaptation and coordinated evolution of plant hydraulic traits. Pablo Sanchez-Martinez, Jordi Martínez-Vilalta, Kyle G. Dexter, Ricardo A. Segovia and Maurizio Mencuccini Affiliations and emails: Pablo Sanchez-Martinez1,2: Email: p.sanchez@creaf.uab.cat Affiliations: (1) CREAF,Cerdanyola del Valles, 08193 Barcelona, Spain; (2) Universitat Autònoma de Barcelona, Cerdanyola del Valles, 08193 Barcelona, Spain Jordi Martínez-Vilalta1,2: Email: Jordi.Martinez.Vilalta@uab.cat Affiliations: (1) CREAF, Cerdanyola del Valles, 08193 Barcelona, Spain (2) Universitat Autònoma de Barcelona, Cerdanyola del Valles, 08193 Barcelona, Spain Kyle G. Dexter3,4: Email: Kyle.Dexter@ed.ac.uk Affiliations: (3) School of GeoSciences, University of Edinburgh, Edinburgh, United Kingdom (4) Royal Botanic Garden Edinburgh, Edinburgh, United Kingdom Ricardo A. Segovia3,5: Email: segovia@ug.uchile.cl Affiliations: (3) School of GeoSciences, University of Edinburgh, Edinburgh, United Kingdom (5) Instituto de Ecología y Biodiversidad (www.ieb chile.cl), Santiago, Chile Maurizio Mencuccini1,6: Email: m.mencuccini@creaf.uab.cat Affiliations: (1) CREAF, Cerdanyola del Valles, 08193 Barcelona, Spain (6) ICREA, Pg. Lluís Companys 23, 08010, Barcelona, Spain Running title: hydraulic traits evolution Keywords: global, plant hydraulics, evolution, phylogenetic comparative methods, phylogenetic mixed models, adaptation, phylogenetic conservatism, evolutionary correlation. 1 Type of Article: Letters Number of words in the abstract: 150 Number of words in the main text: 5145 Number of references: 70 Number of figures: 5 Number of tables: 1 Corresponding Author: Pablo Sanchez-Martinez: p.sanchez@creaf.uab.cat, Carrer Sant Martí de Porres 1-3, 2-4, 08032 Barcelona, Spain. Telephone: +34 644376756. ORCID ID: https://orcid.org/0000-0002-0157-7800. Author contributions: PSM, MM and JMV designed the study, KGD and RSC provided the phylogeny, JMV and MM provided the hydraulics database, PSM analysed data, with input from MM, JMV and KGD, and wrote the first draft of the manuscript. All authors contributed substantially to revisions. Data accessibility statement: if accepted, data will be deposited in a public repository. 2 1 Abstract 2 Hydraulic properties control plant responses to climate and are likely to be under strong selective 3 pressure, but their macro-evolutionary history remains poorly characterized. To fill this gap, we 4 compiled a global dataset of hydraulic traits describing xylem conductivity (Ks), xylem resistance 5 to embolism (P50), sapwood allocation relative to leaf area (Hv) and drought exposure (ψmin), and 6 matched it with global seed plant phylogenies. Individually, these traits present medium to high 7 levels of phylogenetic signal, partly related to environmental selective pressures shaping lineage 8 evolution. Most of these traits evolved independently of each other, being co-selected by the same 9 environmental pressures. However, the evolutionary correlations between P50 and ψmin and 10 between Ks and Hv show signs of deeper evolutionary integration because of functional, 11 developmental or genetic constraints, conforming to evolutionary modules. We do not detect 12 evolutionary integration between conductivity and resistance to embolism, rejecting a hardwired 13 trade-off for this pair of traits. 3 14 Introduction 15 Water transport in plants occurs under negative pressure and is driven by the process of 16 transpiration at the leaf-atmosphere interface, which generates a water potential gradient 17 throughout the plant (cohesion-tension theory) (Dixon 1914). A key source of vulnerability for the 18 water transport system is the formation of xylem embolism, resulting from the breakage of the 19 water columns caused by cavitation (the phase change from liquid water to gas), which reduces 20 hydraulic conductivity and may lead to plant death through hydraulic failure (Tyree & 21 Zimmermann 2002). This process is more likely to occur during drought events, as low water 22 availability results in low soil plant water potentials, and becomes more pronounced also with high 23 temperatures, which provoke an increase in atmospheric evaporative demand (Venturas et al. 24 2017). A wealth of research over the last decades has established that hydraulic failure is a 25 principal mechanism triggering tree mortality under drought (Adams et al. 2017). Therefore, 26 drought and high temperatures, together with other important sources of selection such as freezing 27 temperatures (Zanne et al. 2014), have been considered among the primary forces shaping plant 28 evolution by acting directly on hydraulic traits (Maherali et al. 2004). However, global patterns in 29 the evolution of hydraulic traits remain only partly characterized and their relationship with 30 relevant environmental selective pressures poorly identified. 31 Species differ greatly in their exposure to low water potentials and in their capacity to operate 32 under such conditions. The actual hydraulic risk is normally represented by the hydraulic safety 33 margin (HSM) (Choat et al. 2012). HSM integrates both drought stress exposure at the tissue level, 34 measured as the minimum leaf water potential registered for a given species (ψmin), and resistance 35 to embolism, quantified as the water potential that causes a 50% reduction in stem hydraulic 36 conductivity (P50) (HSM = ψmin - P50). Plants with low (or even negative) safety margins 4 37 experience large amounts of embolism (Choat et al. 2012, 2018). ψmin emerges from the balance 38 between soil water availability, the rate of water loss, and the capacity of the plant transport system 39 to supply water to leaves, and it is thus determined by plant functional properties such as rooting 40 strategy, leaf phenology and stomatal control as well as by abiotic factors such as soil water 41 availability and atmospheric evaporative demand (Bhaskar & Ackerly 2006). Meanwhile, P50 is 42 primarily explained by xylem anatomical features (Venturas et al. 2017). P50 and ψmin are known 43 to co-vary, leading to relatively invariant HSMs at the global scale and to respond to similar 44 environmental selective pressures related to water availability (Choat et al. 2012). For instance, 45 stem P50 has been reported to be negatively related with precipitation for 10 conifer species from 46 different habitats (Brodribb & Hill 1999) and for the gymnosperm genus Callitris (Larter et al. 47 2017) and ψmin has been positively related to variables determining water availability (Bhaskar & 48 Ackerly 2006) and negatively to soil particle size during drought for Great Basin shrubs (Sperry 49 & Hacke 2002). Since the risk of hydraulic failure is likely to be under greater selective pressure 50 than ψmin and P50 per se, these two latter traits are expected to be integrated over the evolutionary 51 history of lineages, specifically meaning that they evolve in a coordinated fashion (i.e., non- 52 independently from each other), representing an evolutionary module. 53 Xylem conductive capacity is another key determinant of hydraulic function, usually quantified as 54 the maximum, stem-specific hydraulic conductivity (Ks). This property has been reported to be 55 positively related to temperature and precipitation at a global scale (He et al. 2020). Because the 56 structural properties of xylem conduits and pit membranes associated with increased embolism 57 resistance (quantified here as P50) are also expected to reduce conductive capacity, a trade-off 58 between P50 and Ks has long been hypothesized (often referred to as the hydraulic safety- 59 efficiency trade-off) (Tyree & Zimmermann 2002). According to this hypothesis, evolutionary 5 60 processes associated with frequent drought occurrence would have driven an increase of xylem 61 resistance to embolism, allowing taxa to bear lower water potentials and maintain water transport 62 at the expense of xylem conductive capacity. In contrast, increased xylem conductivity could have 63 evolved in wetter and warmer environments, where higher water transport was adaptive and 64 selective pressures favouring expensive safety features were weaker (Maherali et al. 2004). 65 Although this trade-off has been shown to be relatively weak across species (e.g. Maherali et al. 66 2004; Gleason et al. 2016), it remains unknown whether it reflects independent responses of each 67 trait to similar selective pressures related to climate conditions and soil properties, or a deeper 68 evolutionary integration. 69 The role of hydraulic conductivity is more nuanced when considered at the whole plant level, 70 where transport capacity needs to match water demand, which is in turn strongly influenced by 71 leaf area (Mencuccini et al. 2019b). Consequently, xylem conductive capacity is frequently 72 expressed in a relativized manner as a measure of hydraulic sufficiency (leaf-specific hydraulic 73 conductivity; Kl, Kl = Ks * Hv, see below) (Tyree & Zimmermann 2002). From this perspective, 74 plants may adapt to drought stress prioritizing supply over demand by reducing the ratio of leaf 75 area relative to cross-sectional sapwood area (i.e., increasing its inverse, the Huber value; Hv) and 76 thus ensuring the maintenance of hydraulic sufficiency under water scarcity. Contrarily, lineages 77 not exposed to drought stress and with no restriction to evolve towards a more conductive xylem 78 may be able to supply water to a higher leaf area by using a relatively low sapwood area, potentially 79 allowing for higher productivity (Mencuccini et al. 2019b). Therefore, we would also expect 80 xylem conductivity and sapwood-to-leaf allocation to be integrated over evolutionary timescales, 81 evolving in a coordinated manner to maintain hydraulic sufficiency (Reich et al. 2003). 6 82 In this study, we aim to elucidate the global macro-evolutionary patterns of hydraulic traits, 83 disentangling (1) the degree to which trait values are evolutionarily conserved along the 84 phylogeny, (2) the extent to which trait conservatism is related to environmentally driven selection 85 and (3) whether traits evolve in a correlated manner because of their independent responses to 86 similar environmental conditions or because of a deeper integration, in which case they may 87 represent evolutionary modules (i.e., a set of traits that co-evolve). We hypothesize that closely 88 related species will have more similar trait values than expected by chance (Losos 2008) and that 89 phylogenetic conservatism will be partly explained by environmental selection (Fig. 1). In 90 addition, we hypothesize that some pairs of traits will show signs of a direct evolutionary 91 relationship (evolutionary modules) reflecting a deep functional, developmental or genetic 92 integration. Specifically, we expect to find three evolutionary modules consistent with previously 93 hypothesized trait coordinations (namely, P50/ψmin, P50/Ks, Ks/Hv). 94 Materials and methods 95 Data sources 96 We extracted detailed hydraulic trait data from a database covering 2027 species (1888 97 angiosperms and 139 gymnosperms), representing 817 genera from 161 families. Most of the data 98 come from a previously published database (Mencuccini et al. 2019b), which was supplemented 99 with the database reported by Liu et al. (2019). Species names were matched against accepted 100 names in The Plant List using the “taxonstand” R package (Cayuela et al. 2012). Then, the 101 “taxonlookup” R package (Pennell et al. 2016) was used to complete species information at the 102 genus, family, order and major evolutionary affiliation (angiosperms vs. gymnosperms) levels. 103 The database covers all major biomes (Fig. S1 in Supporting Information). 104 We used data of four hydraulic traits that were represented across sufficiently large numbers of 7 105 species (N > 550): (1) maximum stem-specific hydraulic conductivity (Ks, kg m-1 MPa-1 s-1) as a 106 measure of xylem conductive capacity; (2) stem water potential at 50% loss of hydraulic 107 conductivity measured in terminal branches (P50, MPa) as a measure of xylem resistance to 108 embolism; (3) branch-based Huber value (Hv; cm2 m-2), defined as the sapwood cross-sectional 109 area to leaf area ratio, as a measure of allocation; and (4) minimum midday leaf water potential 110 recorded for species (ψmin, MPa) as a measure of exposure to drought stress at the tissue level. We 111 also included two additional variables integrating two pairs of the four selected traits, specifically, 112 (5) maximum leaf-specific hydraulic conductivity (Kl, kg m-1 MPa-1 s-1) as the hydraulic capacity 113 per unit leaf area (Ks * Hv) and (6) the hydraulic safety margin (HSM, HSM = ψmin – P50) (Table 114 S1). When multiple measures for one species were available, mean values were used for all traits, 115 except for ψmin, where the absolute minimum was used (c.f. Choat et al. 2012). For all variables, 116 we excluded data from seedlings and studies in greenhouses or experimental gardens, data 117 obtained on roots and leaves (Liu et al. 2019; Mencuccini et al. 2019b) and P50 values 118 corresponding to extreme, r-shaped vulnerability curves, following the same criteria as in Choat 119 et al. (2012). 120 We note that all study traits are subject to methodological uncertainty in their determination and 121 in aggregation to species level, and sample sizes differ among species. Estimates of species- 122 specific ψmin in particular are sample-size dependent and likely biased to an unknown extent for 123 some species. It is likely that the sampling period will miss droughts with a long return interval at 124 some sites. It is also likely that long-lived species (e.g. several gymnosperms) will encounter more 125 severe drought events throughout their lives with consequently greater biases. HSM combines 126 uncertainties in both P50 and ψmin determination, which is problematic because of direct 127 methodological issues in the case of P50 (Jansen et al. 2015) and because of the inherent difficulty 8 128 in characterizing extreme values in the case of ψmin (Head et al. 2012). Finally, in the case of Ks, 129 although it is normalized by sapwood area, it might still depend upon stem size to some degree. 130 Sixteen environmental variables were compiled (11 related to climate and five to soil properties) 131 (Table S1). Climatic variables were extracted from WorldClim (Fick & Hijmans 2017) 132 (www.worldclim.org; accessed on February 2019) except for Moisture Index, which was extracted 133 from the global aridity and potential evapotranspiration (PET) database (Trabucco & Zomer 2019) 134 (http://www.cgiar-csi.org, data accessed on February 2019) at a resolution of 30 arcsec. Soil data 135 were extracted from SoilGrids (Hengl et al. 2017) (https://soilgrids.org/, accessed on February 136 2019) at the same resolution. Occurrences for all species were obtained from the Global 137 Biodiversity Information Facility (https://www.gbif.org/es/ accessed on February 2019) and the 138 Atlas of Living Australia (http://www.ala.org.au. accessed on February 2019) using the “rgbif” 139 (Chamberlain et al. 2020) and the “ALA4R” (Raymond, Vanderwal & Belbin 2014) R packages, 140 respectively. Potentially incorrect species occurrence records where filtered using the 141 “CoordinateCleaner” R package (Zizka et al. 2019). 142 Phylogeny 143 We used a genus-level phylogeny instead of a species-level one to avoid issues with species 144 misidentifications, which are particularly common in the tropics (Baker et al. 2017), and from 145 where a considerable amount of our hydraulics data come. The genera in the phylogeny covered a 146 greater number of species present in our database than the best-sampled species-level phylogeny 147 available (Smith & Brown 2018). Some models, however, were also fitted using the species-level 148 phylogeny from Smith & Brown (2018), to assess the robustness of our results to the taxonomic 149 resolution of our phylogenetic data. To construct the genus-level phylogeny, sequences of the rbcL 150 and matK plastid gene for 707 angiosperm tree genera were obtained from Genbank 9 151 (www.ncbi.nlm.nih.gov/genbank/) building on previous efforts (Dexter & Chave 2016; Neves et 152 al. 2020; Segovia et al. 2020). Sequences were aligned using the MAFFT software (Katoh & 153 Standley 2013). “Ragged ends” of sequences that were missing data for most genera were 154 manually deleted from the alignment. The two chloroplast markers were concatenated, and a 155 maximum likelihood phylogeny for the genera was estimated in the RAxML v8.0.0 software 156 (Stamatakis et al. 2008), on the CIPRES web server (www.phylo.org), using General Time 157 Reversible (GTR) + categorical Gamma (G) model of sequence evolution. The tree was 158 constrained following order-level relationships proposed by the angiosperm Phylogeny Group IV 159 (Chase et al. 2016). Sequences of Nymphaea alba (Nymphaeaceae) were included as an outgroup. 160 The resulting maximum likelihood phylogeny for angiosperms was temporally calibrated using 161 the software treePL (Smith & O’Meara 2012). Age constraints for internal nodes were 162 implemented for most families and orders (Magallón et al. 2015). The rate smoothing parameter 163 (lambda) was set to 10 based on a cross-validation procedure. Finally, the newly-derived 164 angiosperm phylogeny was fused with an existing gymnosperm phylogeny (Leslie et al. 2018). 165 We manually added the genera Gnetum and Ginkgo according to ages found in the literature, 174 166 Ma for the Gnetales (Ran et al. 2018) and 265.2 Ma for Ginkgoaceae (Tank et al. 2015). 167 Statistical analyses 168 All analyses were carried out in R (3.6.0) (R Core team 2019). Some variables were transformed 169 to achieve normality (using absolute values in the case of P50 and ψmin) (Table S1). A Principal 170 Components Analysis (PCA) on the 16 variables was performed using the R package “stats” (R 171 Core team 2019) to reduce the number of axes summarizing environmental variation. The first 172 principal component (PC1) explained 51% of the variance in the environmental data, representing 173 variation in water availability and some related variables such as soil pH, soil clay content, soil 10 174 water content and temperature seasonality, with high values characterizing more humid locations 175 with leached acidic soils characteristic of non-seasonal wet tropical habitats. The second principal 176 component (PC2) explained 20% of the variance, representing variation in energy input, with high 177 values characterizing low solar irradiation, low maximum temperatures and low atmospheric water 178 demand. Finally, the third principal component accounted for 9% of the variance and largely 179 reflected soil depth and, to a lower extent, wind velocity, with high values indicating deeper soils 180 with low sand content and low maximum wind velocities (Table S2, Fig. S2 and Fig. S3). The 181 remaining components explained a low proportion of variance (<7%), so the first three axes were 182 used to characterize the environmental niches of species in the following analyses. 183 Uni-response and bi-response Bayesian phylogenetic mixed models, alternatively including or 184 excluding fixed effects of environmental principal components, major evolutionary affiliation 185 (angiosperm vs. gymnosperm) and their interactions were fitted using the “MCMCglmm” R 186 package (Hadfield 2010) (see Table S3 for models description). All models accounted for the 187 occurrence of multiple measurements in each genus by the inclusion of genus identity as a random 188 effect. Moreover, genus-level phylogenetic relationships were taken into account as a second 189 random effect using the previously presented phylogeny. The inclusion of these random effects 190 allowed us to partition the residual variance from models into three components: the inter-generic 191 variance caused by phylogenetic relationships; the non-phylogenetic, inter-generic variance; and 192 the intra-generic variance. The inter-generic phylogenetic variance quantifies the variability 193 explained by the relationships among taxa as given by our phylogenetic hypothesis and, when 194 divided by the total variance, gives a measure of the phylogenetic signal (λ) (Lynch 1991). Non- 195 phylogenetic inter-generic variance (γ) accounts for the proportion of among-genus variability not 196 explained by the phylogeny, and the intra-generic variance (ρ) provides a measure of the 11 197 proportion of variability caused by intra-generic trait variation (plus any residual error) (Hadfield 198 & Nakagawa 2010) (see Appendix S1 for a more formal description). 199 To partition variances of phylogenetic and non-phylogenetic components, we implemented uni- 200 response models without fixed effects for the six selected hydraulic traits and for the three 201 environmental PCA axes as response variables (Table 1, Table S4 to see non-phylogenetic model 202 variance partitions). To identify relationships between hydraulic traits and environmental PCA 203 axes, we then ran uni-response models with hydraulic traits as response variables and single 204 environmental principal components as fixed effects, both accounting and not accounting for 205 phylogenetic relationships affecting the response trait. To examine the effect of the major split 206 between angiosperms and gymnosperms, additional models included a binary variable describing 207 major evolutionary affiliation and the interaction between affiliation and environment, allowing 208 us to detect statistical differences between angiosperms and gymnosperms in the overall mean 209 values of traits and in their relationships with environmental axes. For each group of nested 210 models, the best fitting one in terms of DIC (Deviance Information Criterion) was selected (Table 211 S5 to see DIC values). Models within 4 DIC units of each other were considered equivalent in 212 terms of fit, and the simplest one was selected. 213 Subsequently, to characterize phylogenetic covariation between the hydraulic traits and between 214 each hydraulic variable and the three environmental principal components, bi-response models 215 were used. In these models, two response variables and their phylogenetic structure were 216 considered simultaneously, resulting in a variance-covariance matrix from which the evolutionary 217 correlation between the two variables could be calculated (Appendix S1). Evolutionary 218 correlations were calculated for all combinations of trait pairs, including and excluding the three 219 environmental components, evolutionary affiliation and their interactions as fixed effects (Fig. S1 12 220 shows data coverage for each combination of traits). Finally, we also estimated evolutionary 221 correlations between traits and single environmental principal components including and 222 excluding evolutionary affiliation as a fixed effect (Table S6 to see all correlations). Bi-response 223 models were also implemented using the species-level phylogeny reported by Smith & Brown 224 (2018) and available in the R package “v.PhyloMaker” (Jin & Qian 2019), to ensure consistency 225 with genus-level results (see Apendix S2). As data availability for the species-level phylogeny was 226 lower, we replicated the bi-response genus-level models using the same reduced dataset to ensure 227 that potential differences between results were not due to different species coverage. We also 228 performed analyses using the species-level phylogeny pruned at the genus-level, to ensure that 229 potential differences between results were not explained by differences in the topologies of our 230 custom-made genus-level phylogeny and the available species-level phylogeny (Appendix S3). 231 Models were specified to achieve convergence while minimizing correlation between iterations 232 (Appendix S1). Marginal variance explained (R2m, variance explained by the fixed effects) and 233 conditional variance explained (R2c, variance explained by both fixed and random effects) were 234 calculated for the uni-response models (Nakagawa & Schielzeth 2013). P-values were calculated 235 for evolutionary correlations following Makowski et al. (2019). 236 Finally, reconstructions of the six traits and the three environmental principal components 237 evolution under a Brownian motion model were mapped along the phylogeny using maximum 238 likelihood ancestral state reconstructions (Schluter et al. 1997) by means of the “Phytools” R 239 package (Revell 2013). 13 240 Results 241 Phylogenetic and non-phylogenetic variances 242 All six selected traits showed a significant phylogenetic signal. The proportion of variance that 243 was explained by the inter-generic phylogenetic structure (λ) ranged from 0.432 (Kl) to 0.745 244 (ψmin) (Table 1). This means that 43.2-74.5% of trait variances can be attributed to relatively deep 245 evolutionary differences among genera, with the rest being attributed to non-phylogenetically 246 correlated inter-generic (γ) and intra-generic (ρ) variances. Intra-generic variances (ρ) ranged from 247 0.189 (ψmin) to 0.532 (Kl), being the second most important variance component in all cases except 248 Kl (where it was the most important), indicating that trait diversification within genera is a 249 substantial generator of global trait variability. Analyses using the species-level phylogeny 250 confirmed that variation within genera also had strong phylogenetic patterns (Appendix S2). 251 Finally, inter-generic, non-phylogenetically related variances (γ) ranged from 0.036 (Kl) to 0.225 252 (P50) and accounted for the lowest proportion of the variance in all cases (Table 1). Phylogenetic 253 mapping of hydraulic traits qualitatively confirmed the findings reported above, showing more 254 gradual changes in highly conserved traits such as ψmin and changes more concentrated at the tips 255 of the phylogeny for variables showing a lower phylogenetic signal, such as Hv, which also 256 showed higher intra-generic variance (Fig. 2, Fig. 3 and Fig. S4). 257 Importantly, the phylogenetic signal of the three environmental principal components was also 258 very high, particularly for PC1, representing water availability (0.820) and PC3, mainly 259 represented by soil depth (0.841) (Table 1, Fig. S5). 260 Environmental drivers of hydraulic trait variability 261 In models that accounted for phylogenetic structure, all hydraulic traits showed significant 262 relationships with at least one of the three environmental axes defined by the PCA (Fig. 4). 14 263 Conditional explained variances (R2c) were notably higher than marginal explained variances 264 (R2m), indicating that accounting for the phylogenetic relationships was crucial to improve model 265 fit (Fig. 4). Consistent with the fact that environmental axes were highly phylogenetically 266 conserved, we also found that the phylogenetic signal of traits (Table 1) diminished when 267 accounting for environmental effects (Fig. 4, lambdas (λ)), thus indicating that environmental 268 conditions explain part of the phylogenetic variance. 269 Xylem resistance to embolism (|P50|) was only negatively related to PC1 (water availability). 270 Minimum water potential at midday (|ψmin|) was negatively related to PC1 and PC2 (declining 271 energy input) and positively to PC3 (soil depth). However, the relationship with PC1 was only 272 significant for angiosperms. Xylem conductivity (Ks) was found to be positively related to PC1 273 and PC3, being negatively related with PC2 only in non-phylogenetic models. Sapwood to leaf 274 area ratio (Hv) was negatively related to PC1 and PC3. The hydraulic safety margin (HSM) was 275 positively related to PC1 and PC2 only for angiosperms and the relationship between HSM and 276 PC3 was only significant (and negative) for non-phylogenetic models. Finally, Leaf-specific 277 conductivity (Kl) was only related to PC2 (negatively) in phylogenetic models, but a positive 278 relationship with PC3 was also found when using non-phylogenetic models (Fig. 4). 279 Evolutionary correlations 280 Significant evolutionary correlations were reported between |ψmin| and |P50| (positive), Ks and Hv 281 (negative), Hv and |P50| (positive) (Fig. 5). These evolutionary correlations were confirmed when 282 the species-level phylogeny was used, which also showed a significant evolutionary correlations 283 between |P50| and Ks (negative), |ψmin| and Hv (positive) and Ks and |ψmin| (negative) (Fig. S6). The 284 emergence of these evolutionary correlations was not explained by the lower number of species 285 available for the species-level phylogeny compared to the genus-level one, nor by differences in 15 286 topology between phylogenies (Appendix S3), so it is likely due to the large amount of 287 phylogenetic covariance between traits within genera. Only two evolutionary correlations between 288 traits remained once environmental factors and major evolutionary affiliation of species were 289 accounted for, coinciding with two of the three hypothesized evolutionary modules. These were 290 the ones involving |P50| and |ψmin| (positive correlation) and Ks and Hv (negative correlation) (Fig. 291 5, Fig. S6). While |P50| and |ψmin| presented a highly conserved covariation pattern, the 292 evolutionary integration between Ks and Hv was less strong. The latter was marginally significant 293 when using the genus-level phylogeny (Fig. 5), but clearly significant when intra-generic 294 phylogenetic covariation between traits was additionally considered by using the species-level 295 phylogeny (Fig, S6). 296 Consistent evolutionary correlations were also observed between certain hydraulic traits and 297 environmental principal components in the bi-response models: Ks was positively correlated with 298 PC1 (water availability), and PC3 (soil depth) while its relationship with PC2 (energy input) was 299 only significant at the genus-level and disappeared when considering major evolutionary 300 affiliation. Hv was negatively correlated with PC1 and PC3; and both |P50| and |ψmin| were 301 negatively correlated with PC1 (Fig. 5). 302 Discussion 303 Conservatism and adaptation in hydraulic trait evolution 304 We found a clear pattern of phylogenetic conservatism for hydraulic traits, suggesting that the 305 legacy of traits found in species’ evolutionary ancestors is an important determinant of trait values 306 in extant species. While we cannot formally rule out Brownian motion evolution operating over 307 long evolutionary timescales as the source of present-day trait variability on the basis of our single 308 trait variance partitioning (Revell et al. 2008), our finding of evolutionary correlations of traits 16 309 with environmental variables indicates a key role for non-random evolutionary processes. 310 Moreover, environmental components explained part of traits’ phylogenetic variance when 311 accounted for as fixed effects (Fig. 4). Therefore, our analyses indicate that adaptive processes 312 have driven the diversification of hydraulic traits, but the prevalent pattern of phylogenetic niche 313 conservatism suggests that evolutionary constraints have limited the range of trait values within 314 lineages. Thus, lineages have been largely tracking environments similar to those their ancestors 315 were already adapted to, retaining ancestral traits because of stabilizing selection (Ackerly 2009), 316 while occasionally adapting to novel environmental conditions. 317 We do observe a wide range of trait values across lineages (including among genera), indicating 318 that they adaptively diverged in deep evolutionary time (Fig. 2, Fig. 3, Fig. S4 and Fig. S5). 319 Further, substantial trait variation can also arise over shorter evolutionary timescales (i.e., in recent 320 evolutionary time) via species adapting to present-day selective pressures, as supported by the 321 significant degree of trait variance at the intra-generic level (Table 1), which also appears to have 322 a phylogenetic component (Fig. S6, Appendix 2). As a result, lineages that have been evolving in 323 dry habitats have adapted to a higher exposure to drought stress by increasing their xylem 324 resistance to embolism, being able to maintain water transport at low water potentials (Choat et al. 325 2012). These species are also selected to ensure water supply to leaves by using a relatively high 326 sapwood area with a low hydraulic conductivity (Mencuccini et al. 2019b). As water become less 327 limiting, lineages are less exposed to low water potentials and are not selected to increase xylem 328 resistance to embolism, while switching their allocation priority to a high leaf area maintained by 329 a smaller area of highly conductive sapwood (Fig. 4). 330 However, substantial variability in species exposure to drought stress within a given environment 331 reflects the fact that plant characteristics such as stomatal control (Brodribb & McAdam 2017), 17 332 deciduousness (Wolfe et al. 2016) or rooting depth (Canadell et al. 1996) are also determining 333 hydraulic trait evolution. This may explain the lack of a relationship between PC1 (water 334 availability) and ψmin and HSM in gymnosperms, a clade mainly represented by Pinaceae and 335 Cupressaceae (Fig. S7) that are known to adopt contrasting strategies under drought. While 336 Pinaceae avoid exposure to low water potentials by closing their stomata and possibly 337 disconnecting from the soil (Poyatos et al. 2018), Cupressaceae tolerate them by presenting a high 338 resistance to embolism (Brodribb et al. 2014). Differences between angiosperms and 339 gymnosperms could also be due to an underestimation of drought stress exposure for long-lived 340 gymnosperms, especially in the case of the highly stress-resistant Cupressaceae, for which the 341 observation window may not have been long enough to adequately capture ψmin. Therefore, 342 different evolutionary processes may be dominant depending on the taxon studied. For instance, 343 xylem resistance has been reported to be extremely labile for the genus Callitris (Larter et al. 344 2017) and to be conserved for Juniperus (Willson et al. 2008), while showing a high canalization 345 for Pinus species (Lamy et al. 2014). It is also worth noting that our global eco-evolutionary 346 overview may be limited by the availability of hydraulic data and its methodological uncertainties, 347 as well as by the difficulty of upscaling traits at the whole-plant level, which remains a challenge 348 (Mencuccini et al. 2019a). 349 Evolutionary modules in hydraulic traits 350 Traits can evolve in an apparently coordinated fashion because of their response to similar selective 351 pressures, but direct relationships between them may also arise from functional, developmental or 352 genetic constraints, conforming to evolutionary modules. We found two sets of traits for which an 353 evolutionary correlation cannot be explained by similar, albeit independent responses to 354 environmental conditions or by fundamental differences between angiosperms and gymnosperms. 18 355 These sets of traits represent a deeper evolutionary integration, confirming two of the three 356 hypothesized evolutionary modules. The first evolutionary module involves species exposure to 357 drought and resistance to embolism (P50/ψmin), and it is strongly conserved over evolutionary 358 scales. The second one involves xylem conductivity and sapwood to leaf area allocation (Ks/Hv), 359 the integration of which appears stronger when quantified in more recent evolutionary time (c.f. 360 results for genus- vs. species-level phylogenies in Fig. 5 and Fig. S6). The third evolutionary 361 module we hypothesized (Ks/P50) appears to be explained exclusively by separate trait responses 362 to similar selective pressures, confirming previous results (Maherali et al. 2004). Therefore, a 363 direct evolutionary trade-off between Ks and P50 can be rejected based on available data, further 364 indicating that Ks and P50 cannot be determined by a single common anatomical feature (e.g., the 365 size distribution of pores in the inter-conduit pit membranes) (Baas et al. 2004). We suggest that 366 Ks and P50 depend on several anatomical properties that may be coordinated under strong selective 367 pressures, but do not necessarily co-evolve when pressures are relaxed over evolutionary 368 timescales. Our results likely reflect the fact that some species may present strategies that do not 369 rely on maximizing xylem conductivity or resistance to embolism, especially when water is not 370 the most limiting resource and survival does not depend on fast resource use (Reich 2014). 371 However, the detailed structural and physiological conditions allowing the independent evolution 372 of these two traits remains to be elucidated. 373 Traits involved in the same evolutionary module are likely to be functionally, developmentally 374 and genetically integrated. Deep functional integrations over evolutionary times can be explained 375 by the need to optimize HSM and Kl under a given environmental context, as the maintenance of 376 positive safety margins and a sufficient hydraulic supply to leaves are likely to be closely linked 377 to survival (Choat et al. 2018) and under a strong stabilizing selection. Therefore, events of 19 378 coordinated directional selection on the involved trait pairs described above might take place over 379 evolutionary times in order to maintain HSM and Kl values close to the adaptive peak. The 380 conservative nature of the relationship between Hv and Ks also reflects broader strategies of 381 convergent evolution integrating hydraulic with photosynthetic and nutrient-use traits as a function 382 of water availability (Hao et al. 2011; Liu et al. 2015). 383 Integration might also be influenced by phylogenetic conservatism in underlying physiological 384 processes and anatomical features. For example, conservatism in stomatal control (Brodribb & 385 McAdam 2017) and leaf phenology (Davies et al. 2013) might explain the evolutionary covariation 386 between ψmin and P50 beyond environmental forcing, with some lineages being able to avoid low 387 water potentials by rapid stomatal closure (Martin-StPaul et al. 2017) or drought-deciduousness 388 (Kolb & Davis 1994). 389 Finally, these functional and developmental integrations may be underpinned by genetic 390 integration, specifically meaning that processes such as genetic correlation (Etterson & Shaw 391 2001), linkage disequilibrium and pleiotropy (Cheverud 1996) might be affecting the anatomical 392 and structural determinants of the traits involved, leading to the observed evolutionary integration. 393 As a result, the evolution of traits in the same module might be genetically constrained (Wagner 394 1996). Further work on the causes and consequences of the evolutionary integration of hydraulic 395 traits, and the meaning of their conservatism through evolutionary time, will be crucial to 396 characterize global trait syndromes and assess species adaptive potential under changing 397 environmental conditions. 20 398 Conclusion 399 Hydraulic traits are under strong selective control and appear to be largely determined by deep- 400 time evolutionary changes driven by adaptation to divergent environmental conditions, in turn 401 limited by evolutionary constraints. We have found evidence for evolutionary integrations not 402 explained by common environmental drivers, conforming to two evolutionary modules: the xylem 403 resistance-exposure module (ψmin/P50), which is highly conserved over evolutionary scales, and 404 the conductivity-allocation module (Hv/Ks), which is more evident in recent evolutionary time. 405 Our results do not support the hypothesized resistance-conductivity module (Ks/P50). The 406 underlying mechanisms shaping these evolutionary modules and their role in species functional 407 and evolutionary diversification remain to be elucidated. More phylogenetically explicit studies of 408 individual clades (including intraspecific genetic, anatomical and functional variation) under 409 different environmental contexts will allow further characterization of plant trait syndromes as a 410 network of integrated units that respond to natural selection. 411 Acknowledgments 412 This work was supported by the Spanish government via competitive grants FUN2FUN 413 (CGL2013-46808-R) and DRESS (CGL2017-89149-C2-1-R).We would like to thank Lucy 414 Rowland, Toby Pennington, Paulo Bittencourt and Mario Marcos do Espirito Santo, and Jarrod 415 Hadfield in particular, for fruitful discussions and insights on previous versions of the manuscript, 416 as well as to Brendan Choat, Steven Jansen and Hui Liu for their previous work on data 417 compilation. P.S.-M. acknowledges an FPU predoctoral fellowship from the Spanish Ministry of 418 Science, Innovation and Universities (grant FPU18/04945). J.M.-V. benefited from an ICREA 419 Academia award. R.A.S. is supported by a Newton International Fellowship from The Royal 420 Society, by Conicyt PFCHA/Postdoctorado Becas Chile/2017 N° 3140189 and by CONICYT PIA 21 421 APOYO CCTE AFB170008. The data derived from the hydraulics database are partly an outcome 422 from a working group funded by the Australian Research Council (ARC) through the Australia- 423 New Zealand Research Network for Vegetation Function. 424 Literature cited 425 Ackerly, D. (2009). Conservatism and diversification of plant functional traits: Evolutionary rates 426 versus phylogenetic signal. Proc. Natl. Acad. Sci., 106, 19699–19706. 427 Adams, H.D., Zeppel, M.J.B., Anderegg, W.R.L., Hartmann, H., Landhäusser, S.M., Tissue, D.T., 428 et al. (2017). A multi-species synthesis of physiological mechanisms in drought-induced tree 429 mortality. Nat. Ecol. Evol., 1, 1285–1291. 430 Baas, P., Ewers, F.W., Davis, S.D. & Wheeler, E.A. (2004). Evolution of xylem physiology. In: 431 (The evolution of plant physiology) {Hemsley, A.R., Poole, L.}. Elsevier, London, 273–295. 432 Baker, T. R., Pennington, R. T., Dexter, K. G., Fine, P. V. A., Fortune-Hopkins, H., Honorio, E. 433 N., et al. (2017). Maximising Synergy among Tropical Plant Systematists, Ecologists, and 434 Evolutionary Biologists. Trends Ecol Evol., 32(4), 258–267. 435 Bhaskar, R. & Ackerly, D.D. (2006). Ecological relevance of minimum seasonal water potentials. 436 Physiol. Plant., 127, 353–359. 437 Brodribb, T. & Hill, R.S. (1999). The importance of xylem constraints in the distribution of conifer 438 species. New Phytol., 143, 365–372. 439 Brodribb, T.J. & McAdam, S.A.M. (2017). Evolution of the stomatal regulation of plant water 440 content. Plant Physiol., 174, 639–649. 22 441 Brodribb, T.J., McAdam, S.A.M., Jordan, G.J. & Martins, S.C.V. (2014). Conifer species adapt to 442 low-rainfall climates by following one of two divergent pathways. Proc. Natl. Acad. Sci., 111, 443 14489–14493. 444 Canadell, J., Jackson, R.B., Ehleringer, J.R., Mooney, H.A., Sala, O.E. & Schulze, E.D. (1996). 445 Maximum rooting depth of vegetation types at the global scale. Oecologia, 108, 583–595. 446 Cayuela, L., Granzow-de la Cerda, Í., Albuquerque, F.S. & Golicher, D.J. (2012). Taxonstand: An 447 r package for species names standardisation in vegetation databases. Methods Ecol. Evol., 3, 1078– 448 1083. 449 Chamberlain S., Barve V., Mcglinn D., Oldoni D., Desmet P., Geffert L. & Ram K. (2020). rgbif: 450 Interface to the Global Biodiversity Information Facility API. R package version 2.2.0. 451 Chase, M.W., Christenhusz, M.J.M., Fay, M.F., Byng, J.W., Judd, W.S., Soltis, D.E., et al. (2016). 452 An update of the Angiosperm Phylogeny Group classification for the orders and families of 453 flowering plants: APG IV. Bot. J. Linn. Soc., 181, 1–20. 454 Cheverud, J.M. (1996). Developmental integration and the evolution of pleiotropy. Am. Zool., 36, 455 44–50. 456 Choat, B., Brodribb, T., Brodersen, C., Duursma, R., López, R. & Medlyn, B. (2018). Triggers of 457 tree mortality under drought drought and forest mortality. Nature, 558, 531–539. 458 Choat, B., Jansen, S., Brodribb, T.J., Cochard, H., Delzon, S., Bhaskar, R., et al. (2012). Global 459 convergence in the vulnerability of forests to drought. Nature, 491, 752–755. 460 Davies, T.J., Wolkovich, E.M., Kraft, N.J.B., Salamin, N., Allen, J.M., Ault, T.R., et al. (2013). 461 Phylogenetic conservatism in plant phenology. J. Ecol., 101, 1520–1530. 23 462 Dexter, K. & Chave, J. (2016). Evolutionary patterns of range size, abundance and species richness 463 in Amazonian angiosperm trees. PeerJ, 2016, 1–14. 464 Dixon, H.H. (1914). Transpiration and the ascent of sap in plants. MacMillan, London, UK. 465 Etterson, J.R. & Shaw, R.G. (2001). Constraint to adaptive evolution in response to global 466 warming. Science, 294, 151–154. 467 Fick, S.E. & Hijmans, R.J. (2017). WorldClim 2: new 1-km spatial resolution climate surfaces for 468 global land areas. Int. J. Climatol., 37, 4302–4315. 469 Gleason, S.M., Westoby, M., Jansen, S., Choat, B., Hacke, U.G., Pratt, R.B., et al. (2016). Weak 470 trad eoff between xylem safety and xylem-specific hydraulic efficiency across the world’s woody 471 plant species. New Phytol., 209, 123–136. 472 Hadfield, J.D. (2010). MCMCglmm for R. J. Stat. Softw., 33. 473 Hadfield, J.D. & Nakagawa, S. (2010). General quantitative genetic methods for comparative 474 biology: Phylogenies, taxonomies and multi-trait models for continuous and categorical 475 characters. J. Evol. Biol., 23, 494–508. 476 Hao, G.Y., Goldstein, G., Sack, L., Holbrook, N.M., Liu, Z.H., Wang, A.Y., et al. (2011). Ecology 477 of hemiepiphytism in fig species is based on evolutionary correlation of hydraulics and carbon 478 economy. Ecology, 92, 2117–2130. 479 He, P., Gleason, S.M., Wright, I.J., Weng, E., Liu, H., Zhu, S., et al. (2020). Growing-season 480 temperature and precipitation are independent drivers of global variation in xylem hydraulic 481 conductivity. Glob. Chang. Biol., 26, 1833–1841. 24 482 Head, A.W., Hardin, J.S. & Adolph, S.C. (2012). Methods for estimating peak physiological 483 performance and correlating performance measures. Environ. Ecol. Stat., 19, 127–137. 484 Hengl, T., Mendes de Jesus, J., Heuvelink, G. B. M., Ruiperez Gonzalez, M., Kilibarda, M., 485 Blagotić, A., et al. (2017). SoilGrids250m: Global gridded soil information based on machine 486 learning. PloS One, 12(2), e0169748. 487 Jansen, S., Schuldt, B. & Choat, B. (2015). Current controversies and challenges in applying plant 488 hydraulic techniques. New Phytol., 205, 961–964. 489 Jin, Y. & Qian, H. (2019). V.PhyloMaker: an R package that can generate very large phylogenies 490 for vascular plants. Ecography, 42, 1353–1359. 491 Katoh, K. & Standley, D.M. (2013). MAFFT multiple sequence alignment software version 7: 492 Improvements in performance and usability. Mol. Biol. Evol., 30, 772–780. 493 Kolb, K.J. & Davis, S.D. (1994). Drought Tolerance and Xylem Embolism in Co-Occurring 494 Species of Coastal Sage and Chaparral. Ecology, 75, 648–659. 495 Lamy, J.B., Delzon, S., Bouche, P.S., Alia, R., Vendramin, G.G., Cochard, H., et al. (2014). 496 Limited genetic variability and phenotypic plasticity detected for cavitation resistance in a 497 Mediterranean pine. New Phytol., 201, 874–886. 498 Larter, M., Pfautsch, S., Domec, J.C., Trueba, S., Nagalingum, N. & Delzon, S. (2017). Aridity 499 drove the evolution of extreme embolism resistance and the radiation of conifer genus Callitris. 500 New Phytol., 215, 97–112. 25 501 Leslie, A.B., Beaulieu, J., Holman, G., Campbell, C.S., Mei, W., Raubeson, L.R., et al. (2018). An 502 overview of extant conifer evolution from the perspective of the fossil record. Am. J. Bot., 105, 503 1531–1544. 504 Liu, H., Gleason, S.M., Hao, G., Hua, L., He, P., Goldstein, G., et al. (2019). Hydraulic traits are 505 coordinated with maximum plant height at the global scale. Sci. Adv., 5, eaav1332. 506 Liu, H., Xu, Q., He, P., Santiago, L.S., Yang, K. & Ye, Q. (2015). Strong phylogenetic signals and 507 phylogenetic niche conservatism in ecophysiological traits across divergent lineages of 508 Magnoliaceae. Sci. Rep., 5, 1–12. 509 Losos, J.B. (2008). Phylogenetic niche conservatism, phylogenetic signal and the relationship 510 between phylogenetic relatedness and ecological similarity among species. Ecol. Lett., 11, 995– 511 1003. 512 Lynch, M. (1991). Methods for the analysis of comparative data in evolutionary biology. 513 Evolution, 45, 1065–1079. 514 Magallón, S., Gómez‐Acevedo, S., Sánchez‐Reyes, L.L. & Hernández‐Hernández, T. (2015). A 515 metacalibrated time‐tree documents the early rise of flowering plant phylogenetic diversity. New 516 Phytol., 207, 437–453. 517 Maherali, H., Pockman, W.T. & Jackson, R.B. (2004). Adaptive variation in the vulnerability of 518 woody plants to xylem cavitation. Ecology, 85, 2184–2199. 519 Makowski, D., Ben-Shachar, M. & Lüdecke, D. (2019). bayestestR: Describing Effects and their 520 Uncertainty, Existence and Significance within the Bayesian Framework. J. Open Source Softw., 521 4, 1541. 26 522 Martin-StPaul, N., Delzon, S. & Cochard, H. (2017). Plant resistance to drought depends on timely 523 stomatal closure. Ecol. Lett., 20, 1437–1447. 524 Mencuccini, M., Manzoni, S. & Christoffersen, B. (2019a). Modelling water fluxes in plants: from 525 tissues to biosphere. New Phytol., 222, 1207-1222. 526 Mencuccini, M., Rosas, T., Rowland, L., Choat, B., Cornelissen, H., Jansen, S., et al. (2019b). 527 Leaf economics and plant hydraulics drive leaf : wood area ratios. New Phytol. 224, 1544-1556. 528 Nakagawa, S. & Schielzeth, H. (2013). A general and simple method for obtaining R2 from 529 generalized linear mixed-effects models. Methods Ecol. Evol., 4, 133–142. 530 Neves, D.M., Dexter, K.G., Baker, T.R., Coelho de Souza, F., Oliveira-Filho, A.T., Queiroz, L.P., 531 et al. (2020). Evolutionary diversity in tropical tree communities peaks at intermediate 532 precipitation. Sci. Rep., 10, 1188. 533 Pennell, M.W., FitzJohn, R.G. & Cornwell, W.K. (2016). A simple approach for maximizing the 534 overlap of phylogenetic and comparative data. Methods Ecol. Evol., 7, 751–758. 535 Poyatos, R., Aguadé, D. & Martínez-Vilalta, J. (2018). Below-ground hydraulic constraints during 536 drought-induced decline in Scots pine. Ann. For. Sci., 75. 537 Ran, J.H., Shen, T.T., Wang, M.M. & Wang, X.Q. (2018). Phylogenomics resolves the deep 538 phylogeny of seed plants and indicates partial convergent or homoplastic evolution between 539 Gnetales and angiosperms. Proc. R. Soc. B Biol. Sci., 285. 540 Raymond, B., Vanderwal, J., & Belbin, L., (2014). ALA4R version 1.01. Atlas of Living 541 Australia. 27 542 R Core Team (2017). R: A language and environment for statistical computing. R Foundation for 543 Statistical Computing, Vienna, Austria. 544 Reich, P. B., Wright, I. J., Cavender‐Bares, J., Craine, J. M., Oleksyn, J., Westoby, M., et al. 545 (2003). The Evolution of Plant Functional Variation: Traits, Spectra, and Strategies. Int. J. Plant 546 Sci., 164(S3), S143–S164. 547 Reich, P.B. (2014). The world-wide “fast-slow” plant economics spectrum: A traits manifesto. J. 548 Ecol., 102, 275–301. 549 Revell, L.J. (2013). Two new graphical methods for mapping trait evolution on phylogenies. 550 Methods Ecol. Evol., 4, 754–759. 551 Revell, L.J., Harmon, L.J. & Collar, D.C. (2008). Phylogenetic signal, evolutionary process, and 552 rate. Syst. Biol., 57, 591–601. 553 Schluter, D., Price, T., Mooers, A. Ø., & Ludwig, D. (1997). Likelihood of ancestor states in 554 adaptive radiation. Evolution, 51(6), 1699–1711. 555 Segovia, R.A., Pennington, R.T., Baker, T.R., De Souza, F.C., Neves, D.M., Davis, C.C., et al. 556 (2020). Freezing and water availability structure the evolutionary diversity of trees across the 557 Americas. Sci. Adv., 6(19), eaaz5373. 558 Smith, S.A. & Brown, J.W. (2018). Constructing a broadly inclusive seed plant phylogeny. Am. J. 559 Bot., 105, 302–314. 560 Smith, S.A. & O’Meara, B.C. (2012). TreePL: Divergence time estimation using penalized 561 likelihood for large phylogenies. Bioinformatics, 28, 2689–2690. 28 562 Sperry, J.S. & Hacke, U.G. (2002). Desert shrub water relations with respect to soil characteristics 563 and plant functional type. Funct. Ecol., 16, 367–378. 564 Stamatakis, A., Hoover, P. & Rougemont, J. (2008). A rapid bootstrap algorithm for the RAxML 565 web servers. Syst. Biol., 57, 758–771. 566 Tank, D.C., Eastman, J.M., Pennell, M.W., Soltis, P.S., Soltis, D.E., Hinchliff, C.E., et al. (2015). 567 Nested radiations and the pulse of angiosperm diversification: increased diversification rates often 568 follow whole genome duplications. New Phytol., 207, 454–467. 569 Trabucco, A. & Zomer, R.. (2019). Global Aridity Index and Potential Evapotranspiration (ET0) 570 Climate Database v2. 571 Tyree M.T. & Zimmermann M.H. (2002). Xylem Structure and the Ascent of Sap. Springer, New 572 York. 573 Venturas, M.D., Sperry, J.S. & Hacke, U.G. (2017). Plant xylem hydraulics: What we understand, 574 current research, and future challenges. J. Integr. Plant Biol., 59, 356–389. 575 Wagner, G.P. (1996). Homologues, natural kinds and the evolution of modularity. Am. Zool., 36, 576 36–43. 577 Willson, C.J., Manos, P.S. & Jackson, R.B. (2008). Hydraulic traits are influenced by phylogenetic 578 history in the drought-resistant, invasive genus Juniperus (Cupressaceae). Am. J. Bot., 95, 299– 579 314. 580 Wolfe, B.T., Sperry, J.S. & Kursar, T.A. (2016). Does leaf shedding protect stems from cavitation 581 during seasonal droughts? A test of the hydraulic fuse hypothesis. New Phytol., 212, 1007–1018. 29 582 Zanne, A.E., Tank, D.C., Cornwell, W.K., Eastman, J.M., Smith, S.A., Fitzjohn, R.G., et al. 583 (2014). Three keys to the radiation of angiosperms into freezing environments. Nature, 506, 89– 584 92. 585 Zizka, A., Silvestro, D., Andermann, T., Azevedo, J., Duarte Ritter, C., Edler, D., et al. (2019). 586 CoordinateCleaner: Standardized cleaning of occurrence records from biological collection 587 databases. Methods Ecol. Evol., 10, 744–751. 30 588 Tables 589 Table 1. Variance partitioning for six hydraulic traits and three environmental principal 590 components related to water availability (PC1), energy input (PC2) and soil depth (PC3). Legend: 591 N: number of species used in each case (for which both phylogenetic and hydraulic data were 592 available), phylogenetic variance (phylogenetic signal, λ), non-phylogenetic inter-generic variance 593 (γ) and intra-generic variance plus measurement error (ρ). Mean and lower and upper 95% credible 594 intervals (HDP) are shown for each component. variable N λ Lower HPD Upper HPD γ Lower HDP Upper HDP ρ Lower HDP Upper HDP Log(|P50|) 868 0.484 0.305 0.697 0.225 0.085 0.360 0.291 0.205 0.368 Log(|ψmin|) 541 0.745 0.572 0.874 0.066 0.000 0.179 0.189 0.129 0.273 log(Ks) 1026 0.515 0.363 0.680 0.086 0.000 0.174 0.399 0.303 0.493 Log(Hv) 1271 0.446 0.291 0.594 0.191 0.097 0.294 0.363 0.276 0.449 HSM 326 0.449 0.201 0.722 0.163 0.000 0.339 0.388 0.246 0.546 Log(Kl) 827 0.432 0.244 0.592 0.036 0.000 0.113 0.532 0.399 0.675 PC1 1911 0.820 0.767 0.870 0.063 0.030 0.099 0.117 0.093 0.139 PC2 1911 0.686 0.599 0.766 0.028 0.000 0.069 0.286 0.230 0.341 PC3 1911 0.841 0.798 0.876 0.007 0.000 0.027 0.152 0.124 0.182 595 31 596 Figure captions 597 Figure 1. Hypotheses and theoretical framework. Double-headed arrows represent potential 598 evolutionary correlation involving key hydraulic traits (xylem conductivity (Ks), xylem resistance 599 to embolism (P50), sapwood allocation relative to leaf area (Huber value, Hv) and drought 600 exposure (ψmin)). HSM refer to Hydraulic Safety Margin, which is the relationship between ψmin 601 and P50, and Kl refer to the hydraulic sufficiency, which is the relationship between Ks and Hv. 602 Lines represent evolutionary relationships tested between pairs of traits. Blue lines represent 603 hypothetical positive relationships between traits, and red lines hypothetical negative ones. Black 604 curved arrows represent traits phylogenetic variance (phylogenetic signal). Each hypothesized 605 coordination between traits is also encircled using long dashed lines and labelled accordingly. 606 Figure 2. Phylogenetic reconstruction of drought exposure (ψmin) and embolism resistance (P50) 607 under a Brownian motion model of evolution. Reconstructions are made log-transformed absolute 608 values in both cases. Families with more than one genus are presented and some of the most 609 important families are highlighted in bold. Gymnosperm families are displayed in grey and 610 angiosperm families in black. 611 Figure 3. Phylogenetic reconstruction of hydraulic conductivity (Ks) and sapwood allocation 612 relative to leaf area (Huber value, Hv). Reconstructions are made on log-transformed data in both 613 cases. Families with more than one genus are presented and some of the most important families 614 are highlighted in bold. Gymnosperm families are displayed in grey and angiosperm families in 615 black. 616 Figure 4. Trait-environment relationships. Relationships between environmental principal 617 components (PC1, which is related to water availability; PC2, which is related to energy input and 32 618 PC3, which is mainly related to soil depth) and hydraulic traits (log-transformed absolute values) 619 accounting for the phylogenetic structure of the hydraulic traits. The best model for each case is 620 displayed, showing the Spermatophyte level relationship (black) or the angiosperm and 621 gymnosperm relationships (red and blue, respectively) when statistically different. Grey dashed 622 lines represent the regression line at the Spermatophyte level without accounting for the 623 phylogenetic structure. Statistically significant (p < 0.05) regression slopes are displayed in bold 624 following the same colour code. Signif. codes: ‘***’: P < 0.001; ‘**’: P < 0.01; ‘*’: P < 0.05 ‘.’: 625 P < 0.1 ‘ ’: P > 0.1. Residual phylogenetic signal (λ) once environmental effects are accounted for 626 in each case is reported when relationships are significant. R2m is the variance explained by the 627 fixed effects and R2c by the fixed and random effects for the phylogenetic mixed models. 628 Figure 5. Evolutionary correlations between hydraulic traits and between traits and environmental 629 principal components (PC1, which is related to water availability; PC2, which is related to energy 630 input and PC3, which is mainly related to soil depth). Environmental variables represent 631 orthogonal PC axes and as such are not correlated. Lines represent significant evolutionary 632 correlations (i.e., when the credible interval for the estimated correlation does not include zero), 633 with the thickness of the line proportional to the strength of the correlation coefficient (also given 634 on the same line). Light red lines represent negative relationships, dark blue lines indicate positive 635 relationships. Significant correlation coefficients between traits when excluding environmental 636 components and evolutionary affiliation as fixed effects are shown in italics, and significant 637 correlation coefficients between traits including environmental components and evolutionary 638 affiliation as fixed effects are shown in bold (in the case of the relationships between 639 environmental axes and traits, only evolutionary affiliation was considered). Dashed lines 640 represent evolutionary correlations that became non-significant when environmental effects and 33 641 major evolutionary affiliations were considered. Pie charts represent phylogenetic signal (dark), 642 inter-generic (medium) and intra-generic (light) variances reported in Table 1, calculated using the 643 maximum number of observations for each case. P-values are also displayed for each coefficient. 644 Signif. codes: ‘***’: P < 0.001; ‘**’: P < 0.01; ‘*’: P < 0.05 ‘.’: P < 0.1 ‘ ’: P > 0.1. 34 645 Figure 1 646 35 647 Figure 2 648 36 649 Figure 3 37 650 Figure 4 651 38 652 Figure 5 653 39 654 Supporting information 655 Figure S1. Whittaker diagrams. 656 Whittaker diagrams for all observations available (once matched with the phylogeny) and 657 observations used for each one of the evolutionary correlations calculation (which has been 658 restricted to those species with complete observations for the two traits and with genus-level 659 phylogenetic information available). 660 40 661 Figure S2. Geographic distribution of the three main environmental principal components. 662 Species-mean coordinates are plotted for each species coloured by their environmental principal 663 components mean values. Thus, some coordinates fall into the sea (presumably species present in 664 both the Palearctic and the Nearctic realms). However, note that environmental variables were 665 calculated for each occurrence of each species separately and then averaged to the species level. 666 PC1 (a) is mainly related to water availability variables, PC2 (b) to energy input and PC3 (b) to 667 soil depth (see Table S2 for a more concrete list of variables and their contribution to each of the 668 three principal component). 41 669 42 670 Figure S3. PCA biplot environment-hydraulic relationships. 671 PCA biplots showing the contributions of the 10 most important environmental variables to the 672 first two principal components, PC1 and PC2 (a) and to PC1 and PC3 (b), colouring species as 673 angiosperms (red circles) or gymnosperms (light blue triangles). Environmental variance 674 explained for each principal component is shown in percentage. log(TS): temperature seasonality 675 (log. transformed); pH: soil pH measured at 60 cm; VPDmax: maximum vapour pressure deficit; 676 Tmax: mean of the monthly maximum temperatures; MAT: mean annual temperature; Clay: clay 677 content in percentage measured at 60cm, log(Wet P): Precipitation of the wettest month (log. 678 Transformed); AP: annual precipitation; AI: aridity index (which is actually a moisture index); 679 sqrt(DQ P): dry quarter precipitation (square root transformed); Soil depth: absolute depth to 680 bedrock, SWC: soil water content at 200cm, Windmax: mean of the monthly maximum wind 681 velocity. 682 43 683 Figure S4. Phylogenetic reconstruction of Hydraulic Safety Margin (HSM) and leaf-specific 684 hydraulic conductivity (log-transformed, log(Kl)) under a Brownian motion model of 685 evolution. 686 Families with more than one genus are presented and some of the most important families are 687 highlighted in bold. Gymnosperm families are displayed in grey and angiosperm ones in black. 44 688 45 689 Figure S5. Phylogenetic reconstruction of the three environmental principal components 690 under a Brownian motion model of evolution. 691 Families with more than one genus are presented and some of the most important families are 692 highlighted in bold. Gymnosperm families are displayed in grey and angiosperm ones in black. 693 PC1 refer to the first environmental principal component, representing variation in water 694 availability and some related variables such as soil pH, soil clay content, soil water content and 695 temperature seasonality, with high values characterizing more humid locations with leached acidic 696 soils characteristic of tropical habitats. PC2 refer to the second environmental principal 697 component, representing variation in energy input, with high values characterizing low solar 698 irradiation, low maximum temperatures and low atmospheric water demand. PC3 refer to the third 699 environmental principal component, largely reflected by soil depth and, to a lower extent, wind 700 velocity, with high values indicating deeper soils with low sand content and low maximum wind 701 velocities. 46 47 702 Figure S6. Evolutionary correlations when using the species-level phylogeny. 703 Evolutionary correlations between hydraulic traits and between traits and environmental variables 704 (environmental variables represent orthogonal PC axes and as such are not correlated) using a 705 species-level phylogeny. Lines represent significant evolutionary correlations (i.e., when the 706 credible interval for the estimated correlation do not include zero), with the thickness of the line 707 proportional to the strength of the correlation coefficient (also given on the same line). Light red 708 lines represent negative relationships, dark blue one’s indicate positive relationships. Significant 709 correlations coefficients between traits when excluding environmental effects and evolutionary 710 affiliation as fixed effects are shown in italics, and significant correlation coefficients between 711 traits including environmental effects and evolutionary affiliation as fixed effects are shown in 712 bold (in the case of the relationships between environmental axis and traits, only evolutionary 713 affiliation was considered). Dashed line represents evolutionary correlation that became non- 714 significant when environmental effects and major evolutionary affiliations were considered. Pie 715 charts represent phylogenetic (dark) and intraspecific (light) variances reported in Appendix 2 (i.e., 716 calculated using the maximum number of observations for each case). P-values are also displayed 717 for each coefficient. Signif. codes: ‘***’: P < 0.001; ‘**’: P < 0.01; ‘*’: P < 0.05 ‘.’: P < 0.1 ‘ ’: P 718 > 0.1. 48 719 49 720 Figure S7. Gymnosperms observations for the relationships between HSM and ψmin with 721 PC1 and the one between HSM and PC2. 722 Species with HSM and ψmin data available are shown coloured by family. PC1 refers to the 723 environmental principal component mainly explained by water availability, PC2 refers to the 724 principal component mainly explained by decreasing energy input. 725 50 726 Table S1. Number of observations variables abbreviation and transformations. 727 Environmental variable and hydraulic traits nomenclature and number of whole dataset and major 728 evolutionary affiliation observations. In the “Transformation” column data transformations are 729 specified, when implemented. Variable Trasformation Abbrebiation Total observations Angiosperms Gymnosperms Potential at the 50% loss of Logaritmic of P50 894 771 123 conductivity the absolute value Maximum stem-specific Logaritmic Ks 1051 951 100 hydraulic conductivity Leaf-specific hydraulic Logaritmic Kl 845 769 76 conductivity Huber value (sapwood area:leaf Logaritmic Hv 1298 1223 75 area ratio) Minimum water potential Logaritmic of ψmin 553 505 48 recorded the absolute value Hydraulic Safety Margin (ψmin- HSM 336 294 42 P50) Precipitation warmest quarter Square root sqrt(WQ P) 1937 1808 129 Precipitation wettest month Logaritmic log(Wet P) 1937 1808 129 Mean of the monthly maximum Tmax 1937 1808 129 temperature Temperature seasonality Logaritmic log(TS) 1937 1808 129 Annual precipitation AP 1937 1808 129 Precipitation driest quarter Square root sqrt(DQ P) 1937 1808 129 Mean annual temperature MAT 1937 1808 129 Aridity index (which is actually a AI 1937 1808 129 moisture index) Solar radiation srad 1937 1808 129 Mean of the monthly maximum windmax 1937 1808 129 wind velocity Maximum vapour pressure deficit VPDmax 1937 1808 129 Absolute depth to bed rock Soil depth 1937 1808 129 pH measured at 60cm pH 1937 1808 129 Clay content in percentage Clay 1937 1808 129 measured at 60cm Sand content in percentage Sand 1937 1808 129 measured at 60cm Soil water content at 200cm SWC 1937 1808 129 51 730 Table S2. Contribution of environmental variables to the three environmental principal 731 components. 732 The highest contribution is highlighted for each variable. Sqrt(WQ P): Precipitation warmest 733 quarter (square root transformed); log(Wet P): Precipitation wettest month (log. Transformed); 734 Tmax: Mean of the monthly maximum temperature; log(TS): Temperature seasonality (log. 735 Transformed); AP: Annual precipitation; sqrt(DQ P): Precipitation driest quarter (square root 736 Transformed); MAT: Mean annual temperature; AI: Aridity index (which is actually a moisture 737 index); srad: Solar radiation; Windmax: Mean of the monthly maximum wind velocity; VPDmax: 738 Maximum vapour pressure deficit; Soil depth: Absolute depth to bedrock; pH: pH measured at 739 60cm; Clay: Clay content in percentage measured at 60cm; Sand: Sand content in percentage 740 measured at 60cm; SWC: Soil water content at 200cm. Contribution Correlation Variable PC1 PC2 PC3 PC1 PC2 PC3 sqrt(WQ P) 7.200 2.273 0.125 0.766 0.269 0.042 log(Wet P) 10.192 0.121 0.031 0.912 0.062 -0.021 Tmax 5.261 16.947 0.022 0.655 -0.734 -0.018 log(TS) 7.906 2.624 2.283 -0.803 0.289 0.181 AP 11.069 0.502 0.588 0.950 0.126 -0.092 sqrt(DQ P) 6.551 5.256 2.105 0.731 0.409 -0.174 MAT 6.749 12.733 0.041 0.742 -0.636 -0.024 AI 8.932 4.820 0.173 0.853 0.391 -0.050 srad 0.077 20.291 11.520 -0.079 -0.803 -0.406 Windmax 5.767 3.142 17.055 -0.686 0.316 -0.495 VPDmax 1.986 19.722 0.920 -0.402 -0.792 0.115 Soil Depth 0.942 1.405 48.430 0.277 -0.211 0.833 pH 8.843 2.569 0.873 -0.849 -0.286 0.112 Clay 6.636 6.592 2.668 0.736 -0.458 -0.196 Sand 4.496 0.825 10.945 -0.606 -0.162 -0.396 SWC 7.393 0.178 2.220 0.776 0.075 -0.178 52 741 Table S3. Reference table for all the models reported in the main text. 742 All models were implemented with and without accounting for the phylogeny. In the fixed 743 structure column, variables to the right of the “~” symbol are response variables, those to the left 744 are predictors. Abbreviations: “env”(1): individual environmental principal component; env(3): 745 three main environmental principal components; trait: individual hydraulic trait; Affiliation: major 746 evolutionary affiliation (angiosperm or gymnosperm), “1” refer to the intercept. Fixed structure Description Phylogeny Number of response variables Results Ref. used env(1) ~ 1 Genus-level Uni-response Table 1, Fig. 4 (pie charts) Phylogenetic signal trait ~ 1 Genus-level Uni-response Table 1, Fig. 4 (pie charts) trait ~ env(1) Genus- Uni-response Fig. 4, Table level S5 trait ~ env(1) + Affiliation Genus- Uni-response Fig. 4, Table Uni-response environment models level S5 trait ~ env(1) * Affiliation Genus- Uni-response Fig. 4, Table level S5 trait , env(1) ~ 1 Genus- Bi-response Fig. 5, Table level S6 trait , env(1) ~ 1 + Affiliation Genus- Bi-response Fig. 5, Table level S6 trait , trait ~ 1 Genus- Bi-response Fig. 5, Table level S6 Evolutionary correlations trait , trait ~ 1 + Affiliation Genus- Bi-response Fig. 5, Table level S6 trait , trait ~ 1 + env(3) Genus- Bi-response Fig. 5, Table level S6 trait , trait ~ 1 + env(3) + Affiliation Genus- Bi-response Fig. 5, Table level S6 53 trait , trait ~ 1 + env(3) * Affiliation Genus- Bi-response Fig. 5, Table level S6 env(1) ~ 1 Species- Uni-response Appendix 2, level Fig. S6 Phylogenetic signal trait ~ 1 Species- Uni-response Appendix 2, level Fig. S6 trait , env(1) ~ -1 Species- Bi-response Appendix 2, level Fig. S6 trait , env(1) ~ -1 + Affiliation Species- Bi-response Appendix 2, level Fig. S6 Evolutionary correlations trait , trait ~ -1 Species- Bi-response Appendix 2, level Fig. S6 trait , trait ~ -1 + env(3) * Affiliation Species- Bi-response Appendix 2, level Fig. S6 747 54 748 Table S4. Non-phylogenetic model’s variance partition. 749 Mean non-phylogenetic inter-generic (γ) and non-phylogenetic intra-generic (ρ) variance in non- 750 phylogenetic models without fixed effects. Note that phylogenetic variance (λ) is 0, as the 751 phylogenetic effect was not considered. variable Phylogenetic (λ) Inter-generic (γ) Intra-generic (ρ) HSM 0 0.490 0.510 Log(Hv) 0 0.514 0.486 Log(Kl) 0 0.280 0.720 Log(Ks) 0 0.459 0.541 Log(|ψmin|) 0 0.621 0.379 Log(|P50|) 0 0.636 0.364 PC1 0 0.787 0.213 PC2 0 0.483 0.517 PC3 0 0.641 0.359 55 752 Table S5. Uni-response models description. 753 DICs and explained variances for phylogenetic and non-phylogenetic uni-response models. The 754 fixed formula is shown in each case. DICs for the phylogenetic models are shown. “NP” refer to 755 non-phylogenetic models (i.e., only including genus contingency as random effect) explained 756 variances. R2c refer to the conditional and R2m refers to the marginal explained variances. 757 Abbreviations: Ks: Xylem conductivity; P50: xylem resistance to embolism; Hv: sapwood 758 allocation relative to leaf area; ψmin: drought exposure, HSM: hydraulic safety margin; Kl: and 759 sufficiency; PC1: water availability; PC2: energy input and PC3: soil depth; Affiliation: 760 evolutionary affiliation (angiosperm or gymnosperm). Fixed effects formula DIC R2m R2c NP R2m NP R2c HSM ~ 1 1181 0 0.612 0 0.49 HSM ~ PC1 * Affiliation 1133 0.253 0.657 0.301 0.554 HSM ~ PC1 + Affiliation 1138 0.211 0.647 0.268 0.535 HSM ~ PC1 1139 0.065 0.625 0.026 0.519 HSM ~ PC2 * Affiliation 1133 0.246 0.623 0.28 0.509 HSM ~ PC2 + Affiliation 1142 0.184 0.658 0.237 0.54 HSM ~ PC2 1143 0.035 0.619 0.031 0.495 HSM ~ PC3 * Affiliation 1146 0.172 0.624 0.23 0.542 HSM ~ PC3 + Affiliation 1147 0.176 0.634 0.235 0.541 HSM ~ PC3 1150 0.006 0.618 0.02 0.528 log_Hv ~ 1 3147 0 0.641 0 0.514 log_Hv ~ PC1 + Affiliation 3013 0.166 0.549 0.187 0.479 log_Hv ~ PC1 3013 0.155 0.536 0.184 0.477 log_Hv ~ PC1 * Affiliation 3014 0.168 0.551 0.188 0.478 log_Hv ~ PC2 * Affiliation 3058 0.028 0.662 0.014 0.51 log_Hv ~ PC2 3060 0.002 0.646 0.001 0.509 log_Hv ~ PC2 + Affiliation 3060 0.028 0.657 0.013 0.509 log_Hv ~ PC3 + Affiliation 3066 0.045 0.615 0.04 0.477 log_Hv ~ PC3 3066 0.022 0.601 0.032 0.48 log_Hv ~ PC3 * Affiliation 3068 0.045 0.614 0.043 0.475 log_Kl ~ 1 2348 0 0.47 0 0.28 log_Kl ~ PC1 * Affiliation 2250 0.067 0.482 0.077 0.299 log_Kl ~ PC1 2252 0.002 0.447 0.002 0.282 log_Kl ~ PC1 + Affiliation 2254 0.069 0.463 0.077 0.296 56 log_Kl ~ PC2 2250 0.016 0.41 0.033 0.254 log_Kl ~ PC2 * Affiliation 2251 0.089 0.439 0.093 0.283 log_Kl ~ PC2 + Affiliation 2251 0.084 0.434 0.088 0.276 log_Kl ~ PC3 2250 0.002 0.454 0.008 0.297 log_Kl ~ PC3 + Affiliation 2252 0.072 0.462 0.077 0.304 log_Kl ~ PC3 * Affiliation 2253 0.075 0.47 0.077 0.304 log_Ks ~ 1 2795 0 0.608 0 0.459 log_Ks ~ log_Hv 2079 0.116 0.614 0.181 0.494 log_Ks ~ log_Hv * Affiliation 2081 0.166 0.64 0.23 0.506 log_Ks ~ log(|P50|) 1581 0.041 0.634 0.074 0.499 log_Ks ~ log(|P50|) * Affiliation 1583 0.111 0.666 0.118 0.51 log_Ks ~ PC1 2670 0.028 0.603 0.042 0.475 log_Ks ~ PC1 + Affiliation 2670 0.084 0.631 0.089 0.479 log_Ks ~ PC1 * Affiliation 2670 0.091 0.635 0.09 0.479 log_Ks ~ PC2 2694 0.003 0.602 0.01 0.447 log_Ks ~ PC2 + Affiliation 2694 0.057 0.623 0.055 0.453 log_Ks ~ PC2 * Affiliation 2696 0.058 0.625 0.056 0.456 log_Ks ~ PC3 2685 0.01 0.59 0.028 0.462 log_Ks ~ PC3 + Affiliation 2685 0.063 0.618 0.069 0.462 log_Ks ~ PC3 * Affiliation 2687 0.062 0.617 0.07 0.463 log(|ψmin|) ~ 1 828 0 0.812 0 0.621 log(|ψmin|) ~ log(|P50|) 431 0.279 0.738 0.299 0.71 log(|ψmin|) ~ log(|P50|) * Affiliation 432 0.29 0.738 0.314 0.703 log(|ψmin|) ~ PC1 * Affiliation 688 0.228 0.788 0.303 0.686 log(|ψmin|) ~ PC1 692 0.211 0.763 0.302 0.679 log(|ψmin|) ~ PC1 + Affiliation 692 0.229 0.781 0.302 0.68 log(|ψmin|) ~ PC2 * Affiliation 724 0.099 0.854 0.107 0.692 log(|ψmin|) ~ PC2 725 0.055 0.843 0.094 0.691 log(|ψmin|) ~ PC2 + Affiliation 725 0.099 0.854 0.104 0.687 log(|ψmin|) ~ PC3 + Affiliation 793 0.07 0.841 0.047 0.649 log(|ψmin|) ~ PC3 793 0.016 0.827 0.036 0.644 log(|ψmin|) ~ PC3 * Affiliation 795 0.072 0.841 0.051 0.65 log(|P50|) ~ 1 1426 0 0.71 0 0.636 log(|P50|) ~ PC1 * Affiliation 1396 0.193 0.635 0.23 0.605 log(|P50|) ~ PC1 1397 0.069 0.617 0.097 0.588 log(|P50|) ~ PC1 + Affiliation 1397 0.194 0.636 0.231 0.606 log(|P50|) ~ PC2 * Affiliation 1402 0.116 0.725 0.148 0.635 log(|P50|) ~ PC2 + Affiliation 1403 0.108 0.719 0.141 0.631 log(|P50|) ~ PC2 1403 0.001 0.694 0.002 0.623 log(|P50|) ~ PC3 * Affiliation 1394 0.107 0.728 0.141 0.637 log(|P50|) ~ PC3 + Affiliation 1397 0.105 0.725 0.144 0.636 57 log(|P50|) ~ PC3 1397 0.003 0.699 0.002 0.628 PC1 ~ 1 6985 0 0.891 0 0.787 PC2 ~ 1 6723 0 0.697 0 0.483 PC3 ~ 1 4707 0 0.85 0 0.641 761 Table S6. Evolutionary correlations and DICs for bi-response models. 762 Mean of the evolutionary correlation (Cor), credible interval (lower and upper HDP) and p-value 763 reported by bi-response models. The fixed formula is shown in each case. Models are ordered by 764 DIC values (from lower to higher) for each set of nested models (same response variables). 765 Statistically significant evolutionary correlations are highlighted in bold and marginally significant 766 in italics. Abbreviations: Ks: Xylem conductivity; P50: xylem resistance to embolism; Hv: 767 sapwood allocation relative to leaf area; ψmin: drought exposure, HSM: hydraulic safety margin; 768 Kl: and sufficiency; PC1: water availability; PC2: energy input and PC3: soil depth; Affiliation: 769 evolutionary affiliation (angiosperm or gymnosperm). Var. 1 Var. 2 Fixed formula Cor Lower HDP Upper HDP p-value DIC log(Hv) log(|ψmin|) (log(Hv), log(|ψmin|)) ~ 1 + (PC1 + PC2 + PC3) * :Affiliation -0.094 -0.681 0.436 0.752 1403 log(Hv) log(|ψmin|) (log(Hv), log(|ψmin|)) ~ 1 + PC1 + PC2 + PC3 -0.100 -0.636 0.486 0.736 1403 log(Hv) log(|ψmin|) (log(Hv), log(|ψmin|)) ~ 1 + Affiliation + PC1 + PC2 + PC3 -0.117 -0.624 0.462 0.664 1404 log(Hv) log(|ψmin|) (log(Hv), log(|ψmin|)) ~ 1 + Affiliation 0.222 -0.405 0.801 0.494 1571 log(Hv) log(|ψmin|) (log(Hv), log(|ψmin|)) ~ 1 0.217 -0.405 0.798 0.509 1571 log(Hv) PC1 (log(Hv), PC1) ~ 1 -0.796 -0.913 -0.662 0.000 7475 log(Hv) PC1 (log(Hv), PC1) ~ 1 + Affiliation -0.792 -0.910 -0.657 0.000 7475 log(Hv) PC2 (log(Hv), PC2) ~ 1 + Affiliation 0.145 -0.160 0.462 0.397 7296 log(Hv) PC2 (log(Hv), PC2) ~ 1 0.156 -0.141 0.473 0.363 7296 log(Hv) PC3 (log(Hv), PC3) ~ 1 + Affiliation -0.558 -0.747 -0.363 0.000 5971 log(Hv) PC3 (log(Hv), PC3) ~ 1 -0.565 -0.737 -0.367 0.000 5972 log(Ks) log(Hv) (log(Ks)_log(Hv)) ~ 1 + PC1 + PC2 + PC3 -0.423 -0.805 -0.016 0.077 3843 log(Ks) log(Hv) (log(Ks)_log(Hv)) ~ 1 + Affiliation + PC1 + PC2 + PC3 -0.423 -0.795 -0.005 0.079 3843 cbind(log_(Ks)_log(Hv)) ~ 1 + (PC1 + PC2 + PC3) * log(Ks) log(Hv) -0.421 -0.827 -0.012 0.085 3853 Affiliation log(Ks) log(Hv) (log(Ks), log(Hv)) ~ 1 -0.600 -0.868 -0.271 0.008 4058 log(Ks) log(Hv) (log(Ks), log(Hv)) ~ 1 + Affiliation -0.588 -0.879 -0.247 0.010 4059 log(Ks) log(|ψmin|) (log(Ks), log(|ψmin|)) ~ 1 + (PC1 + PC2 + PC3) * Affiliation 0.019 -0.538 0.566 0.934 1571 log(Ks) log(|ψmin|) (log(Ks), log(|ψmin|)) ~ 1 + Affiliation + PC1 + PC2 + PC3 -0.013 -0.535 0.569 0.966 1572 log(Ks) log(|ψmin|) (log(Ks), log(|ψmin|)) ~ 1 + PC1 + PC2 + PC3 0.008 -0.512 0.555 0.970 1572 58 log(Ks) log(|ψmin|) (log(Ks), log(|ψmin|)) ~ 1 -0.080 -0.683 0.487 0.785 1768 log(Ks) log(|ψmin|) (log(Ks), log(|ψmin|)) ~ 1 + Affiliation -0.066 -0.675 0.507 0.823 1768 log(Ks) log(|P50|) (log(Ks), log(|P50|)) ~ 1 + Affiliation + PC1 + PC2 + PC3 -0.046 -0.517 0.399 0.851 2489 log(Ks) log(|P50|) (log(Ks), log(|P50|)) ~ 1 + PC1 + PC2 + PC3 -0.098 -0.510 0.324 0.648 2489 log(Ks) log(|P50|) (log(Ks), log(|P50|)) ~ 1 + (PC1 + PC2 + PC3) * Affiliation -0.019 -0.475 0.408 0.917 2495 log(Ks) log(|P50|) (log(Ks), log(|P50|)) ~ 1 -0.317 -0.665 0.055 0.114 2596 log(Ks) log(|P50|) (log(Ks), log(|P50|)) ~ 1 + Affiliation -0.274 -0.639 0.165 0.211 2596 log(Ks) PC1 (log(Ks), PC1) ~ 1 + Affiliation 0.339 0.093 0.578 0.010 6343 log(Ks) PC1 (log(Ks), PC1) ~ 1 0.351 0.110 0.594 0.009 6343 log(Ks) PC2 (log(Ks), PC2) ~ 1 + Affiliation -0.275 -0.570 0.003 0.069 6315 log(Ks) PC2 (log(Ks), PC2) ~ 1 -0.292 -0.576 -0.007 0.048 6315 log(Ks) PC3 (log(Ks), PC3) ~ 1 + Affiliation 0.412 0.144 0.683 0.008 5216 log(Ks) PC3 (log(Ks), PC3) ~ 1 0.426 0.185 0.706 0.005 5216 log(|ψmin|) PC1 (log(|ψmin|), PC1) ~ 1 + Affiliation -0.776 -0.908 -0.624 0.000 2636 log(|ψmin|) PC1 (log(|ψmin|), PC1) ~ 1 -0.784 -0.907 -0.620 0.000 2636 log(|ψmin|) PC2 (log(|ψmin|), PC2) ~ 1 -0.232 -0.566 0.129 0.214 2744 log(|ψmin|) PC2 (log(|ψmin|), PC2) ~ 1 + Affiliation -0.249 -0.594 0.124 0.195 2745 log(|ψmin|) PC3 (log(|ψmin|), PC3) ~ 1 + Affiliation 0.115 -0.236 0.452 0.539 1916 log(|ψmin|) PC3 (log(|ψmin|), PC3) ~ 1 0.126 -0.208 0.480 0.471 1917 log(|P50|) log(Hv) (log(|P50|), log(Hv)) ~ 1 + PC1 + PC2 + PC3 0.060 -0.423 0.507 0.818 1834 log(|P50|) log(Hv) (log(|P50|), log(Hv)) ~ 1 + Affiliation + PC1 + PC2 + PC3 0.061 -0.413 0.511 0.779 1835 log(|P50|) log(Hv) (log(|P50|), log(Hv)) ~ 1 + (PC1 + PC2 + PC3) * Affiliation 0.053 -0.404 0.531 0.831 1843 log(|P50|) log(Hv) (log(|P50|), log(Hv)) ~ 1 0.414 0.052 0.790 0.070 1927 log(|P50|) log(Hv) (log(|P50|), log(Hv)) ~ 1 + Affiliation 0.385 -0.043 0.777 0.121 1927 (log(|P50|), log(|ψmin|)) ~ 1 + (PC1 + PC2 + PC3) * log(|P50|) log(|ψmin|) 0.571 0.213 0.860 0.015 745 Affiliation log(|P50|) log(|ψmin|) (log(|P50|), log(|ψmin|)) ~ 1 + Affiliation + PC1 + PC2 + PC3 0.552 0.168 0.863 0.030 751 log(|P50|) log(|ψmin|) (log(|P50|), log(|ψmin|)) ~ 1 + PC1 + PC2 + PC3 0.556 0.168 0.833 0.015 751 log(|P50|) log(|ψmin|) (log(|P50|), log(|ψmin|)) ~ 1 + Affiliation 0.702 0.428 0.917 0.000 834 log(|P50|) log(|ψmin|) (log(|P50|), log(|ψmin|)) ~ 1 0.683 0.386 0.914 0.006 834 log(|P50|) PC1 (log(|P50|), PC1) ~ 1 + Affiliation -0.725 -0.885 -0.537 0.000 4476 log(|P50|) PC1 (log(|P50|), PC1) ~ 1 -0.714 -0.875 -0.524 0.000 4477 log(|P50|) PC2 (log(|P50|), PC2) ~ 1 0.050 -0.316 0.460 0.803 4487 log(|P50|) PC2 (log(|P50|), PC2) ~ 1 + Affiliation 0.015 -0.344 0.400 0.978 4488 log(|P50|) PC3 (log(|P50|), PC3) ~ 1 + Affiliation -0.110 -0.456 0.208 0.570 3677 log(|P50|) PC3 (log(|P50|), PC3) ~ 1 -0.097 -0.426 0.229 0.573 3678 770 59 771 Appendix S1. Supplementary methods. 772 Phylogenetic mixed model description 773 Phylogenetic mixed models are commonly used in quantitative genetics (the so called “animal” 774 model), being useful for comparative analyses as they allow to incorporate a range of variance 775 structures for the random effects, including shared ancestry through a phylogeny (Housworth et 776 al. 2004). The general model structure is defined as follows: 777 𝑦 = 𝜇 + 𝛽𝑥 + 𝑝 + 𝑔 + 𝑒 (1) 778 Were μ is the grand mean, interpreted as the root ancestor state, β is the slope for the covariate x 779 (fixed effect, in green), p and g are the variability caused by the genus-level phylogeny and the 780 genus contingency effects (random effects, in red), and e is the residual error (Housworth et al. 781 2004; Villemereuil & Nakagawa 2014). Both fixed (β) and random (r, which is p + g) effects and 782 the residuals (e) are expected to come from a multivariate normal distribution as it follows: 𝛽 𝛽0 𝐵 0 0 783 [ 𝑟 ] ~ 𝑁 ([ 0 ] , [ 0 𝐺 0 ]) (2) 𝑒 0 0 0 𝑅 784 Where β is the fixed effect parameter to estimate, β0 is the prior means for the fixed effects with 785 prior (co)variance matrix B, and G and R are the expected (co)variances of the random effects and 786 the residuals respectively (Hadfield 2010; Hadfield & Nakagawa 2010). G and R are unknown, 787 and must be estimated from the data by assuming they are structured in a way that can be 788 parametrized by few parameters, as it has been exemplified below for the G case: 𝑉𝐺 ⊗ 𝐴𝐺1 0 789 𝐺= [ 1 ] (3) 0 𝑉𝐺2 ⊗ 𝐴𝐺2 790 Were the (co)variance matrices (V) are matrices with one parameter to be estimated per response 60 791 variable and the structured matrices (A) refer to the phylogenetic structure (A G1) and genus 792 contingency (AG2). The Kronecker product (⊗) allows for possible dependence between random 793 effects (Hadfield 2010; Hadfield & Nakagawa 2010). 794 In multi-response models, the (co)variance matrix of the previous equation is reformulated 795 including the covariance estimates in the off-diagonal and the respective variances in the diagonal 796 as follows: 𝜎𝑢2 𝜎𝑢1 ,𝑢2 797 𝑉𝐺1 = [ 1 ] (4) 𝜎𝑢2 ,𝑢1 𝜎𝑢22 798 Where 𝜎2u1 is the variance for the first response variable (V1) and σ2u2 the variance for the second 799 response variable (V2), while σu1,u2 and σu2,u1 are the same covariance estimate (C). 800 Phylogenetic indexes calculation 801 The phylogenetic signal or phylogenetic heritability it is calculated as follows ( Villemereuil & 802 Nakagawa 2014): 𝜎𝑝2 803 𝜆= (5) 𝜎𝑝 + 𝜎𝑔2 + 𝜎𝑒2 2 804 Where 𝜎𝑃2 is the variance of the phylogenetic effect (VG1), 𝜎𝑔2 is the variance of the cross-genus 805 effect (VG2) and 𝜎𝑒2 is the residual error (Villemereuil & Nakagawa 2014). Cross-genera variance 806 (i.e. non-phylogenetic variation among genera or genus lability) has been calculated as follows: 𝜎𝑔2 807 𝛾= (6) 𝜎𝑝 + 𝜎𝑔2 + 𝜎𝑒2 2 808 And finally, intra-genus variability including measurement error has been calculated as follows: 𝜎𝑒2 809 𝜌= (7) 𝜎𝑝 + 𝜎𝑔2 + 𝜎𝑒2 2 810 Note also that γ + ρ + λ = 1 (Housworth et al. 2004). The three indexes were calculated for the 61 811 whole Markov chain random effects and residual samples (once burned and thinned), so the output 812 is a statistical distribution from which the mean and 95% credible intervals can be calculated. 813 Phylogenetic covariation calculation 814 From the phylogenetic variances and covariance in equation 4, the evolutionary correlation 815 between response variables can be calculated as follows (Villemereuil 2012): 𝜎𝑢2 ,𝑢1 816 𝑟𝑒𝑣 = (8) √𝜎𝑢21 · 𝜎𝑢22 817 Model specifications 818 “MCMCglmm“ implements a Bayesian approach, estimating the posterior distribution of 819 parameters, from which 95% credible intervals can be obtained (Hadfield 2010). We set 820 independent normal prior distributions for fixed effects and non-informative Inverse-Gamma prior 821 distributions for random effects and residual variances (Villemereuil & Nakagawa 2014). Less 822 informative expanded priors were also used, and highly similar results were obtained. 823 Uni-response models random effects variance priors were set as V = 1, nu = 0.002. For bi-response 824 models, the random effects variances priors were set as V = diag(2)/2, nu = 2. To achieve 825 convergence, each model was run for 8,000,000 iterations with a 1,000,000 burn-in and a thinning 826 interval of 4,000, reaching an effective sample size between 1,000 and 2,000 in all estimated 827 parameters. When models did not converge, we increased the number of iterations until 828 convergence were achieved. Thinning intervals and the final number of iterations were 829 progressively increased until autocorrelations between samples were found to be <0.1. 830 Convergence of all models was assessed by plots of chain mixing and by the Heidenberg stationary 831 test as a diagnostic. All reported models had a low degree of autocorrelation between iterations 832 and passed the convergence diagnostic, both for fixed and random effects (i.e., the sampled chains 62 833 were stationary). 834 Literature cited 835 Hadfield, J.D. (2010). MCMCglmm for R. J. Stat. Softw., 33. 836 Hadfield, J.D. & Nakagawa, S. (2010). General quantitative genetic methods for comparative 837 biology: Phylogenies, taxonomies and multi-trait models for continuous and categorical 838 characters. J. Evol. Biol., 23, 494–508. 839 Housworth, E.A., Martins, E.P. & Lynch, M. (2004). The Phylogenetic Mixed Model. Am. Nat., 840 163, 84–96. 841 Villemereuil, P. (2012). How to use the MCMCglmm R package. 842 Villemereuil, P. & Nakagawa, S. (2014). General Quantitative Genetic Methods for Comparative 843 Biology. In: (Modern Phylogenetic Comparative Methods and Their Application in 844 Evolutionary Biology: Concepts and Practice), {[Garamszegi, L.Z.]}. Springer, New York, 845 USA. 287-303. 846 Appendix S2. Species-level phylogenetic analyses. 847 Species-level phylogeny was obtained by pruning the phylogenetic tree reported by Smith & 848 Brown (2018) available in the R package “v.PhyloMaker” (Jin & Qian 2019) by using the “ape” 849 R package (Paradis & Schliep 2018) only keeping species with hydraulic data available in each 850 case, obtaining the same number of observations compared to the genus-level analyses. Some bi- 851 response models implemented using the genus-level phylogeny where also conducted using the 852 species-level phylogeny. As we had only one value per specie, no extra random effect was 853 included, so variance partition was reduced to phylogenetic signal calculation. 63 854 Phylogenetic signal results 855 Variance partitioning for the six hydraulic traits and three environmental principal components 856 related to water availability (PC1), energy input (PC2) and soil depth (PC3). Legend: N: number 857 of species used in each case (for which both phylogenetic and hydraulic data were available), 858 phylogenetic variance (phylogenetic signal, λ) and non-phylogenetic intraspecific variance plus 859 measurement error (ρ). Mean and lower and upper 95% credible intervals (HDP) are shown for 860 each component. variable N λ λ lower HPD λ upper HPD ρ ρ lower HDP ρ Upper HDP HSM 195 0.456 0.228 0.680 0.544 0.320 0.772 Log(Hv) 842 0.654 0.539 0.774 0.346 0.226 0.461 Log(Kl) 616 0.610 0.456 0.753 0.390 0.247 0.544 Log(Ks) 763 0.681 0.569 0.792 0.319 0.208 0.431 Log(|ψmin|) 358 0.876 0.799 0.940 0.124 0.060 0.201 log_negP50 693 0.709 0.594 0.817 0.291 0.183 0.406 PC1 1329 0.963 0.951 0.975 0.037 0.025 0.049 PC2 1329 0.845 0.796 0.889 0.155 0.111 0.204 PC3 1329 0.907 0.882 0.934 0.093 0.066 0.118 861 64 862 Evolutionary correlations results 863 Mean of the evolutionary correlation (Cor), credible interval (lower and upper HDP) and p-value 864 reported by bi-response models. Statistically significant evolutionary correlations are highlighted 865 in bold and marginally significant in italics. In the fixed structure column, variables to the right of 866 the “~” symbol are response variables, those to the left are predictors. Abbreviations: “env”(1): 867 individual environmental principal component; env(3): three main environmental principal 868 components; trait: individual hydraulic trait; Affiliation: major evolutionary affiliation 869 (angiosperm or gymnosperm). Fixed structure Var. 1 Var. 2 Cor Lower Upper HDP p-value HDP trait, trait ~ 1 + env(3) * Affiliation log(Hv) log(|ψmin|) 0.134 -0.365 0.695 0.641 trait, trait ~ 1 log(Hv) log(|ψmin|) 0.607 0.261 0.915 0.014 trait, env(1) ~ 1 log(Hv) PC1 -0.807 -0.908 -0.699 0.000 trait, env(1) ~ 1 + Affiliation log(Hv) PC1 -0.816 -0.922 -0.714 0.000 trait, env(1) ~ 1 log(Hv) PC2 -0.090 -0.334 0.191 0.495 trait, env(1) ~ 1 + Affiliation log(Hv) PC2 -0.092 -0.376 0.164 0.501 trait, env(1) ~ 1 + Affiliation log(Hv) PC3 -0.492 -0.689 -0.304 0.000 trait, env(1) ~ 1 log(Hv) PC3 -0.493 -0.691 -0.304 0.000 trait, trait ~ 1 + env(3) * Affiliation log(Ks) log(Hv) -0.630 -0.851 -0.359 0.000 trait, trait ~ 1 log(Ks) log(Hv) -0.589 -0.815 -0.348 0.000 trait, trait ~ -1 + env(3) * Affiliation log(Ks) log(|ψmin|) -0.217 -0.663 0.226 0.349 trait, trait ~ 1 log(Ks) log(|ψmin|) -0.366 -0.703 0.012 0.090 trait, trait ~ -1 + env(3) * Affiliation log(Ks) log(|P50|) -0.236 -0.579 0.172 0.223 trait, trait ~ 1 log(Ks) log(|P50|) -0.420 -0.674 -0.104 0.015 trait, env(1) ~ 1 log(Ks) PC1 0.225 0.000 0.421 0.043 trait, env(1) ~ 1 + Affiliation log(Ks) PC1 0.225 0.006 0.452 0.067 trait, env(1) ~ 1 + Affiliation log(Ks) PC2 -0.185 -0.434 0.065 0.160 trait, env(1) ~ 1 log(Ks) PC2 -0.196 -0.439 0.092 0.155 trait, env(1) ~ 1 + Affiliation log(Ks) PC3 0.106 -0.132 0.350 0.395 trait, env(1) ~ 1 log(Ks) PC3 0.105 -0.147 0.338 0.423 trait, env(1) ~ 1 + Affiliation log(|ψmin|) PC1 -0.734 -0.861 -0.599 0.000 trait, env(1) ~ 1 log(|ψmin|) PC1 -0.743 -0.868 -0.590 0.000 trait, env(1) ~ 1 + Affiliation log(|ψmin|) PC2 -0.266 -0.573 0.041 0.127 trait, env(1) ~ 1 log(|ψmin|) PC2 -0.254 -0.567 0.040 0.118 trait, env(1) ~ 1 + Affiliation log(|ψmin|) PC3 0.215 -0.032 0.462 0.101 65 trait, env(1) ~ 1 log(|ψmin|) PC3 0.223 -0.032 0.453 0.097 trait, trait ~ -1 + env * Affiliation log(|P50|) log(Hv) 0.211 -0.256 0.663 0.429 trait, trait ~ 1 log(|P50|) log(Hv) 0.622 0.370 0.839 0.001 trait, trait ~ -1 + env(3) * Affiliation log(|P50|) log(|ψmin|) 0.773 0.582 0.926 0.000 trait, env(1) ~ 1 log(|P50|) log(|ψmin|) 0.794 0.636 0.923 0.000 trait, env(1) ~ 1 + Affiliation log(|P50|) PC1 -0.466 -0.658 -0.254 0.000 trait, env(1) ~ 1 log(|P50|) PC1 -0.465 -0.661 -0.257 0.000 trait, env(1) ~ 1 + Affiliation log(|P50|) PC2 0.022 -0.250 0.305 0.902 trait, env(1) ~ 1 log(|P50|) PC2 0.032 -0.225 0.343 0.837 trait, env(1) ~ 1 + Affiliation log(|P50|) PC3 -0.147 -0.417 0.102 0.262 trait, env(1) ~ 1 log(|P50|) PC3 -0.144 -0.390 0.118 0.279 870 871 Literature cited 872 Jin, Y., & Qian, H. (2019). V.PhyloMaker: an R package that can generate very large phylogenies 873 for vascular plants. Ecography, 42(8), 1353–1359. 874 Paradis, E., & Schliep, K. (2018). ape 5.0: an environment for modern phylogenetics and 875 evolutionary analyses in R. Bioinformatics 35: 526-528. 876 Smith, S. A., & Brown, J. W. (2018). Constructing a broadly inclusive seed plant phylogeny. Am. 877 J. Bot., 105(3), 302–314. 66 878 Appendix S3. Evolutionary correlations reported by genus-level phylogenetic models using 879 observations available for the species-level phylogeny and evolutionary correlations 880 reported by species-level phylogeny pruned at the genus level. 881 For bivariate models including two traits as response variable, only models without fixed effects 882 and models including the three environmental components and its interaction with major 883 evolutionary affiliation (angiosperm or gymnosperm) were implemented. 884 Significant evolutionary correlations (i.e., when the credible interval for the estimated correlation 885 do not include zero) reported by models using a genus-level phylogeny including only observations 886 available for the species-level phylogenetic analyses to check for effects of the different species 887 coverage between phylogenies. Mean of the evolutionary correlation (Cor), credible interval 888 (lower and upper HDP) and p-value reported by bi-response models. In the fixed structure column, 889 variables to the right of the “~” symbol are response variables, those to the left are predictors. 890 Abbreviations: “env”(1): individual environmental principal component; env(3): three main 891 environmental principal components; trait: individual hydraulic trait; Affiliation: major 892 evolutionary affiliation (angiosperm or gymnosperm). Fixed structure var1 var2 Cor Lower HDP Upper HDP p-value trait, env(1) ~ 1 + Affiliation log(Hv) PC1 -0.779 -0.926 -0.634 0.000 trait, env(1) ~ 1 log(Hv) PC1 -0.787 -0.921 -0.647 0.000 trait, env(1) ~ 1 + Affiliation log(Hv) PC3 -0.499 -0.749 -0.250 0.001 trait, env(1) ~ 1 log(Hv) PC3 -0.510 -0.749 -0.260 0.001 trait, trait ~ 1 log(Ks) log(Hv) -0.501 -0.850 -0.101 0.049 trait, trait ~ 1 + env(3) * Affiliation log(Ks) log(Hv) -0.603 -0.861 -0.297 0.003 trait, trait ~ 1 log(Ks) log(|P50|) -0.394 -0.709 -0.022 0.054 trait, env(1) ~ 1 + Affiliation log(Ks) PC2 -0.316 -0.602 -0.011 0.049 trait, env(1) ~ 1 log(Ks) PC2 -0.341 -0.618 -0.056 0.024 trait, env(1) ~ 1 + Affiliation log(Ks) PC3 0.350 0.045 0.633 0.031 trait, env(1) ~ 1 log(Ks) PC3 0.355 0.065 0.651 0.021 trait, env(1) ~ 1 + Affiliation log(|ψmin|) PC1 -0.779 -0.926 -0.623 0.000 67 trait, env(1) ~ 1 log(|ψmin|) PC1 -0.783 -0.928 -0.621 0.000 trait, trait ~1 log(|P50|) log(Hv) 0.495 0.126 0.816 0.014 trait, trait ~ 1 log(|P50|) log(|ψmin|) 0.485 0.065 0.836 0.054 trait, trait ~ 1 +env(3) * Affiliation log(|P50|) log(|ψmin|) 0.598 0.233 0.888 0.008 trait, env(1) ~ 1 + Affiliation log(|P50|) PC1 -0.628 -0.863 -0.394 0.000 trait, env(1) ~ 1 log(|P50|) PC1 -0.618 -0.831 -0.374 0.001 893 894 Significant evolutionary correlations (i.e., when the credible interval for the estimated correlation 895 do not include zero) reported by models using a species-level phylogeny pruned at genus-level to 896 check for effects of differences in the topology between phylogenies. Mean of the evolutionary 897 correlation (Cor), credible interval (lower and upper HDP) and p-value reported by bi-response 898 models. In the fixed structure column, variables to the right of the “~” symbol are response 899 variables, those to the left are predictors. Abbreviations: “env”(1): individual environmental 900 principal component; env(3): three main environmental principal components; trait: individual 901 hydraulic trait; Affiliation: major evolutionary affiliation (angiosperm or gymnosperm). Fixed structure var1 var2 Cor Lower HDP Upper HDP p-value trait, env(1) ~ 1 + Affiliation log(Hv) PC1 -0.817 -0.924 -0.685 0.000 trait, env(1) ~ 1 log(Hv) PC1 -0.824 -0.936 -0.696 0.000 trait, env(1) ~ 1 + Affiliation log(Hv) PC3 -0.439 -0.705 -0.158 0.008 trait, env(1) ~ 1 log(Hv) PC3 -0.451 -0.711 -0.159 0.012 trait, trait ~ 1 log(Ks) log(Hv) -0.535 -0.844 -0.178 0.018 trait, trait ~ 1 + env(3) * Affiliation log(Ks) log(Hv) -0.626 -0.877 -0.330 0.002 trait, trait ~ 1 log(Ks) log(|P50|) -0.398 -0.749 -0.016 0.069 trait, env(1) ~ 1 + Affiliation log(Ks) PC2 -0.326 -0.626 -0.013 0.046 trait, env(1) ~ 1 log(Ks) PC2 -0.332 -0.688 -0.037 0.068 trait, env(1) ~ 1 log(Ks) PC3 0.334 0.005 0.647 0.045 trait, env(1) ~ 1 + Affiliation log(|ψmin|) PC1 -0.774 -0.924 -0.614 0.000 trait, env(1) ~ 1 log(|ψmin|) PC1 -0.783 -0.933 -0.619 0.000 trait, trait ~1 log(|P50|) log(Hv) 0.505 0.131 0.814 0.024 trait, trait ~ 1 log(|P50|) log(|ψmin|) 0.493 0.066 0.858 0.054 trait, trait ~ 1 +env(3) * Affiliation log(|P50|) log(|ψmin|) 0.591 0.177 0.906 0.034 trait, env(1) ~ 1 + Affiliation log(|P50|) PC1 -0.609 -0.856 -0.349 0.001 trait, env(1) ~ 1 log(|P50|) PC1 -0.595 -0.831 -0.342 0.001 902 68