Xylem conduits have lignified walls to resist crushing pressures. The thicker the double-wall (T) relative to its diameter (D), the greater the implosion safety. Having safer conduits may incur higher costs and reduced flow, while having less resistant xylem may lead to catastrophic collapse under drought. Although recent studies have shown that conduit implosion commonly occurs in leaves, little is known about how leaf xylem scales T versus D to trade off safety, flow efficiency, mechanical support, and cost.

We measured T and D in > 7,000 conduits of 122 species to investigate how T versus D scaling varies across clades, habitats, growth forms, leaf and vein sizes.

As conduits become wider, their double-cell walls become proportionally thinner, resulting in a negative allometry between T and D. That is, narrower conduits, which are usually subjected to more negative pressures, are proportionally safer than wider ones. Higher implosion safety (i.e higher T/D ratios) was found in asterids, arid habitats, shrubs, small leaves, and minor veins.

Despite the strong allometry, implosion safety does not clearly trade off with other measured leaf functions, suggesting that implosion safety at whole-leaf level cannot be easily predicted solely by individual conduits’ anatomy.

Study species

We sampled 122 species from the University of California Botanical Garden at Berkeley (‘UCBG’, 37.87, -122.23; Berkeley, CA, USA). Species selection maximized phylogenetic coverage, and included species with different growth forms, leaf sizes, and habitats. Because this garden collection often has a single or just a few individuals per species, our sampling approach was limited to a few branches (> 1 m long) collected from a single individual (woody species), or to a few leaves collected from 1-5 individuals (herbaceous species). Samples were collected in the morning, re-cut under water, re-hydrated overnight, and then used for the measurement of leaf anatomical and functional traits. As most of those traits are destructive, we used different leaves for each trait. 

  Anatomical measurements

Fresh and mature leaves (3-4 per species) were cut into 1 cm² sections, fixed in FAA (formalin acetic acid), and embedded in paraffin blocks. Sections included different leaf portions (base, middle, apex, petiole) and vein orders. Transverse cross-sections of 8–10 μm of thickness were cut with a microtome (Leica, RM2265), stained using the Johansen’s safranin-O and fast green method, and mounted in permanent glass slides using Cytoseal 60 medium (Richard-Allan Scientific). Sections were observed under a light microscope (Leica, DM 2000), and photographed (at 20 –100× objectives) with a camera control unit (Nikon, DS-Fi1). We selected 6-8 images per species and manually measured the maximum and minimum lumen diameters (D_max and D_min, μm) and the double-cell wall thickness (T, μm) on all or, at most, 10 adjacent conduits per image using ImageJ (https://imagej.nih.gov/). In vascular bundles with more than 10 conduits we systematically selected conduits to cover the range of conduit sizes observed in each picture.

Functional traits

Implosion safety ratio (T/D) was calculated as the ratio between T and D_max  for each conduit. Higher T/D values are assumed to be associated with greater implosion safety. Per this assumption, we estimated a theoretical critical implosion pressure (P_cri, MPa), i.e. the critical pressure above which a cell wall collapses. P_cri should be a function of T/D, but the specific equation depends on what type of stress (radial or hoop) induces the implosion. Because the prevalent type of stress causing leaf conduits to collapse is not known precisely, we used two different mechanical models of P_cri.

The first model (Blackman et al., 2010), is based on Timoshenko's equation for an isolated and perfectly rounded pipe under negative pressure ( and considers that hoop forces are the main underlying stress inducing collapse. In the absence of per-species data, P_cri1 was calculated assuming that the radial elastic modulus of xylem conduits ranges from 100 (P_cri1 low) to 300 MPa (P_cri1 high) (Blackman et al., 2010). To estimate how much the leaf conduits cross-sectional shape deviates from a rounded shape, we calculate the conduit ovality as O = (D_max - D_min)/(D_max + D_min), where D_min is the minimum conduit lumen diameter in μm (Ikeda et al., 2013). An ovality value of 0 indicates a perfectly rounded conduit. Conduits with O < 0.005 still behave as a cylindrical pipe under compressive forces, but above this value conduit shape starts interfering in P_cri1 (Ikeda et al., 2013). Note that O was not directly used to adjust P_cri1.

The second model ( Hacke et al., 2001), considers the double-wall between neighboring conduits as a flat solid plate of finite D_max and effectively infinite length. It also assumes that radial stresses occurring in the common wall between an embolized and a water-filled conduit are the main cause of implosion, while hoop stresses are negligible (Sperry and Hacke, 2004). Cell wall strength is unknown for the species evaluated. In wood xylem, cell wall strength ~ 40-80 MPa, but it could be lower for leaf xylem. Therefore, we calculated P_cri2 assuming that W ranges from 10 (P_cri2 low) to 80 MPa (P_cri2 high). Although both mechanical models described above are approximations and may not give accurate pressures for collapse, they likely account for enough variation across species to provide useful insights.

Flow efficiency was quantified as the maximum leaf hydraulic conductance (Kleaf_max, mmol m^-2 s^-1 MPa^-1), measured on 4-10 leaves per species using the evaporative flux method (Sack and Scoffoni, 2012) with a pressure-drop flow meter (Melcher et al., 2012). In this method, a transpiring leaf was firmly connected to a tube running to a water source, which was then placed in series with a resistance tube (PEEK tubing, VWR, Radnor, PA, USA) of known hydraulic conductance. To accelerate the evaporation process, leaf samples were placed over a box fan and below a light source. Once a steady-state flow rate was achieved, pressures across the resistance tube were recorded, for at least 10 minutes, using two pressure transducers (model PX26-001GV, Omega Engineering, Norwalk, CT, USA) interfaced to a data logging system (U6 USB, Labjack, Lakewood, CO, USA). Next, the leaf was disconnected from the tubing system, and placed in a sealable plastic bag for about 20 minutes for water potential equilibration. Final leaf water potential was measured using a pressure chamber (model 1505D, PMS, Albany, OR, USA); and leaf area (LA, cm²) was obtained using a flatbed scanner and the leafarea macro (https://github.com/bblonder/leafarea) in ImageJ software version 1.53t (https://imagej.nih.gov/). Finally, Kleaf_max  normalized by leaf area and corrected for leaf temperature was calculated following Sack and Scoffoni (2012). Kleaf_max describes how much water flows across the leaf in response to a water potential gradient between the leaf and the surrounding atmosphere, hence higher values indicate higher flow efficiency. LA obtained as above, was used to classify species as microphyllous (LA ≤ 20.25 cm²), mesophyllous (20.25 < LA < 45 cm²), and macrophyllous (LA ≥ 45 cm²).

Mechanical support was quantified as the leaf flexural modulus of elasticity (ε, MN m^-2). To obtain ε, we performed 3-point bending tests using a Universal Testing Machine (Test stand ES30, Mark-10, Copiague, NY, USA) on 3-4 leaves per species. Leaves were placed in the UTM machine with their longitudinal axis parallel to the bending fixture. During the bending test, the force (force gauges M5-5 and M5-20, Mark-10, Copiague, NY, USA) and the displacement (travel display ESM0001, Mitutoyo, Aurora, IL, USA) were recorded, and then used to produce force-displacement plots. After each test, leaf width and thickness at the bending point and the span length between the two bending fixtures were measured with a digital caliper, and then used for the calculation of ε following Read et al. (2005). Higher values of ε (i.e. stiffer leaves), indicate higher mechanical support. In 14 species, leaves were too small to be properly attached to the bending fixture and/or too flexible to produce detectable bending forces, so ε was not measured.

Construction cost was estimated using two different proxies: leaf mass per area (LMA, g m⁻²) and total volume of veins per area (VTotV, mm³ mm⁻²). LMA describes the total amount of resources invested in constructing each unit of leaf area. To obtain LMA, 3-5 leaves per species were scanned to obtain the leaf area and oven-dried at 50 °C for 48 hours to determine their dry mass. LMA was then calculated as leaf dry mass divided by leaf area (Pérez-Harguindeguy et al., 2016). VTotV describes the total volume of veins per unit of leaf area and it is a reasonable proxy of the construction cost of conduits per se. To calculate VTotV, we first obtained leaf cleared images of all species, except for Nymphaea spp. and Aucuba japonica Thunb., because it was impossible to obtain a clear image of their venation networks. Next, we used the LeafVeinCNN app in Matlab to calculate the total volume of veins (Xu et al., 2021) assuming that veins have a cylindrical shape. Finally, we divided the total volume of veins by leaf area to obtain VTotV. Higher LMA and VTotV values reflect higher construction cost.

  Species habitat

We inferred habitat for each species based on its current geographic range. To retrieve the geographic range, we used occurrence data from BIEN (Botanical Information and Ecology Network) and GBIF (Global Biodiversity Information Facility, https://www.gbif.org/) databases. Occurrence records were manually cleaned by removing duplicated, outdated (pre-1950), or suspect geographical coordinates (outside the species natural habitat). Next, the mean annual precipitation (MAP) for each coordinate was extracted from the Worldclim database (https://www.worldclim.org/data/bioclim.html) at 2.5’ resolution (~ 5 km). For each species, we averaged the MAP for all occurrence points and then classified the species habitat as hydric (MAP > 2000 mm), mesic (500 mm < MAP ≤ 2000 mm) or arid (MAP ≤ 500 mm).

Statistical analysis

To investigate how leaf conduits’ T and D_max scale to each other, we log₁₀-transformed both variables and then used the ‘smatr’ R-package version 3 (Warton et al., 2012) to fit standardized major axis (SMA) regression models. We used the functions sma(log(T) ~ log(D), slope.test = 1) to test if the coefficient b for all species together was significantly different from one. To investigate whether the coefficients b (log(T) ~ log(D) * group) and a (log(T) ~ log(D) + group) differed across groups (species, clades, habitats, growth forms, leaf sizes, and vein orders), we used Likelihood ratio (𝜆) and Wald (W) statistics, respectively (Warton et al., 2012). We also performed Kruskal Wallis tests followed by pairwise Wilcox tests with Benjamini and Hochberg p-value adjustment method to test for differences in anatomical traits across those groups.  Currently, the ‘smatr’ package does not support multiple regressions (Warton et al., 2012). Thus, to investigate any potential influence of both leaf sizes and species in the T ~ D scaling relationship across vein orders we fitted a standard mixed model regression using the function “lme” from the ‘nlme’ R-package. Leaf area was model as fixed effect, while vein orders nested within species, and species nested within genus as random effects, i.e. lme(log(T) ~ log(D) * log(LA), random= ~1|genus/species/vein orders). We used the ‘nlme’ function ‘anova()’ to assess the significance of each predictor variable effect and the function ‘residplot()’ from the ‘predictmeans’ R-package to perform residual analysis.

To further examine phylogenetically driven variation in implosion safety, we built a phylogenetic tree for all species using ‘V.PhyloMaker2’ R-package (Jin and Qian, 2022) and the ‘GBOTB.extended.WP.tre’ mega-tree. Next, we used the function ‘fastAnc’from the ‘phytools’ R-package to perform a fast estimation of the ancestral states for a, b, and T/D. We also conducted Blomberg’s K tests (Blomberg et al., 2003) using the ‘phytools’ function ‘phylosig’ to test for phylogenetic signals in those same traits (i.e. the tendency for related species to resemble each other, more than they resemble species drawn at random from a tree). K < 1 indicates faster trait evolution than expected under a Brownian model (weaker phylogenetic signal), while K > 1 indicates slower trait evolution (stronger phylogenetic signal).

To investigate possible trade-offs among leaf functions, we used two complementary approaches. First, we carried out a principal component analysis (PCA) using the ‘prcomp’ function in R. Prior to the PCA, we centered and z-transformed all traits to improve comparability among them and reduce bias towards traits with higher variance. We used the broken stick method for estimating the number of principal components to be retained. Our PCA analysis was carried out with 106 out of the 122 studied species, as we removed 16 species with missing data for VTotV and ε. Second, we ran ordinary least-squares regression models to test for pairwise trade-offs between implosion safety (response variable) and the other leaf traits (predictor variables). We also regressed P_cri1 and P_cri2 values to assess the relationship between our two different models of conduit collapse.

All analyses were carried out using the R version 4.3.1. Code to reproduce analyses is available at https://github.com/ilamatos/xylem_implosion_safety. Regressions or differences were considered significant if P <0.05.