Climate and hydraulic traits interact to set thresholds for liana viability

TRY meta-­analysis

We used the TRY plant trait database²⁷ to identify traits that show systematic differences between the tree and liana growth forms, as a way to narrow the scope of the rest of the analysis. We chose traits to represent major trade­offs within the “economic spectrum” framework, which places plants along a spectrum of strategies from acquisitive, fast return on investment to conservative, slow return on investment according to key functional trait values³⁰. We narrowed traits to those that had observations for at least four tree and liana species. We then compiled our dataset using the following steps during November and December 2019. For each trait, we downloaded the dataset for all species available globally and averaged the observations of the trait to the species level to avoid statistical biases introduced in our growth form comparison due to a high density of observations in a few commercially valuable species. We matched the species ID number with the most frequently used growth form identifier using the TRY “growth form” trait and kept the species with the most frequent identifier of “tree,” “liana,” or “woody vine.” We subsetted the data to keep only species with a majority of observations ascribed to the tree and liana growth forms (i.e., no herbaceous species, ferns, etc.), resulting in observations for 44,222 total species. Finally, we filtered the dataset of 44,222 species by hand to remove species misclassified as trees or lianas; species occurring entirely in temperate to boreal biomes; species from all gymnosperm lineages except the order Gnetales; and entries for taxonomic classifications broader than the genus level (e.g., taxonomic families). We found that hydraulic functional traits in the TRY database (i.e., K_s,max and P₅₀) show systematic differences between growth forms (Supplementary Fig. 1; Supplementary Tables 3 and 4), while there is mixed evidence for differences in the acquisitiveness of trees and lianas in terms of stem anatomical traits (Supplementary Fig. 1; Supplementary Tables 3 and 4) and leaf functional traits (Supplementary Fig. 6; Supplementary Tables 3 and 4), and no evidence of differences between tropical trees and lianas with respect to root functional traits (Supplementary Fig. 7; Supplementary Tables 3 and 4).

Extended meta­-analysis

We conducted an additional literature search to supplement the hydraulic trait observations from the TRY database. The additional literature search served two purposes: (1) to fill a major gap identified during our TRY analysis in terms of liana trait observations, and (2) to address the methodological inconsistency of measuring K_s,max and P₅₀ on liana branches shorter than the longest vessel, which incorrectly measures K_s,max and P₅₀ without accounting for end wall resistivity^59,60.

We conducted a literature search using Web of Science and Google Scholar. We searched the following phrases in combination with “liana:” “hydraulic conductivity,” “hydraulic trait,” “hydraulic efficiency,” and “hydraulic K.” Of the literature we found, we kept only the studies that met the following criteria: (1) reported K_s,max measurements for lianas, (2) measured K_s,max instead of computing K_s,max from xylem conduit dimensions, (3) measured K_s,max on sunlit, terminal branches of mature individuals or saplings, and (4) measured K_s,max on a branch longer than the longest vessel. We considered the authors to have used a branch length longer than maximum vessel length if the authors measured or reported maximum vessel length for the species and a longer branch was used. Because the best methodological practice for measuring P₅₀, especially in species with long vessels, is currently a matter of debate, we additionally removed all observations of P₅₀ > ­0.75. This filtering was performed to reduce the probability that falsely high (i.e., less negative) P₅₀ values were retained in our analysis because of improper measurement technique and is consistent with the P₅₀ filtering performed by Trugman et al.⁶¹. Improper measurement technique is a particular concern for lianas, whose wide and long vessels require cautious implementation of the traditional measurement techniques developed for trees. We note that retaining all liana P₅₀ observations (i.e., not filtering out observations > −0.75) results in a significant difference between trees and lianas (Mann­–Whitney test statistic = 1029, n_tree = 61, n_liana = 46, p < 0.05). However, the effect size remains relatively small, indicating that even when retaining unrealistically high liana P₅₀ values, the difference between liana and tree P₅₀ is ecologically of only moderate significance (Glass’ Δ = 0.47). When possible, we manually inspected vulnerability curves of each species and removed strongly r-­shaped curves, but corresponding hydraulic safety margins were not available for a quantitative determination. We applied the same criteria to the observations in the TRY database, combined the observations from TRY and from our additional literature search, and averaged the observations to the species level. This resulted in a total of 154 species with hydraulic trait observations matching our criteria, of which 51 species were lianas and 103 species were trees.

A list of the sources of our measurements is available in Supplementary Table 5^{35,62,63,64,65,66,67,68,69,70,71,72}.

Statistical analysis

For both the TRY analysis and the extended meta­analysis, we compared the tree and liana growth forms using two methods. First, we used two­-sided Mann–­Whitney U-­tests, which test whether observations between groups are drawn from the same distribution. We used Mann–­Whitney U-tests rather than t-­tests because the distributions of most traits violate the normality assumption of t-­tests. This approach is consistent with a recent pantropical meta-analysis comparing liana and tree functional trait distributions⁷³.

Second, we computed Glass’ ∆, a measure of effect size, which describes the magnitude of the difference between groups compared with the variation within the reference group^74,75. The Glass’ ∆ was chosen rather than Cohen’s d because the standard deviation of each group is substantially different for several traits, including hydraulic traits^74,75. To avoid biasing our interpretation of the statistics by considering only one growth form as the reference group, we computed and present the test statistic and 95% confidence intervals resulting from using both the tree growth form (subscript “T”) and liana growth form (subscript “L”) as the reference group (Supplementary Table 2; Supplementary Table 4). Throughout the text, we present the statistics computed using the tree as the reference group for two reasons. First, we were interested in the degree to which lianas differ from the well-parameterized tree plant functional types in dynamic vegetation models. Second, because lianas are often parameterized using data from early successional tropical trees⁵⁵, we were interested in considering the degree to which the distribution of liana trait values differs from the distribution of tree trait values characterizing the plant functional types in which lianas are traditionally categorized.

All statistical analyses were conducted in the R statistical environment⁷⁶. Mann–­Whitney U-tests were conducted using the “stats” package and Glass’ ∆ statistics were computed using the “effectsize” package⁷⁷.

Competition model

We modified the single-­tree model originally developed by Trugman et al.³⁹ to represent a single liana-­tree pairing. The purpose of the original model developed by Trugman et al. is to calculate annual net primary production (A_net) of a single temperate tree under defined climatic conditions and morphological and physiological parameters, with A_net becoming the input to a subsequent model describing tree drought recovery. Briefly, the model couples water transport using the Shinozaki pipe model⁴¹ and the Ball-Berry model of stomatal conductance⁴² and whole-plant photosynthesis using the Farquhar photosynthesis model⁴⁰. The amount of water moving through the plant depends on soil water availability (soil water potential, Ψ); the hydraulic path length and xylem area of fine roots, stem, and petioles; and the water demand imposed on the tree by the atmosphere (vapor pressure deficit, VPD). Mathematically, this can be written with the following set of equations. First, the flow, F (mmol s⁻¹), throughout a plant element is computed by integrating the hydraulic conductivity per unit of xylem area (K) from one end of the pipe continuum with water potential ψ₁ (MPa) to the other with water potential ψ₂, which can be expressed by the differences in the Kirchhoff transforms as

where a (m²) is the xylem area of the element and L (m) the pipe length. The element conductivity (K, mmol m⁻¹ s⁻¹ MPa⁻¹) decreases as stem water potential falls as a result of embolism. A logistic function of the form

where b1 is the slope of the percent loss of conductivity (PLC) curve and b2 is P₅₀, is used to represent the loss of conductivity as water potential becomes more negative, and thus ϕ (mmol m⁻¹ s⁻¹) is a function of the maximum whole-plant hydraulic conductivity, K_max (mmol m⁻¹ s⁻¹ MPa⁻¹). The assumptions of our pipe model (i.e., constant xylem area, a, with branching and path length, L, that is representative of the whole path from roots to leaves) allows us to approximate an individual tree or liana with an effective element conductivity for the entire path. This is in contrast to stem-specific hydraulic conductivity (K_s,max, mmol m⁻¹ s⁻¹ MPa⁻¹), which is commonly measured in the field on terminal branches and does not account for the tapering of vessel elements in branches. Therefore, K_max is distinct from K_s,max.

If we neglect changes in water storage, F is constant throughout the hydraulic continuum. Then, water flow from the roots to the stem is modeled as

where a_root and a_stem are the cross-sectional xylem area of the root system and the cross-sectional area of the xylem, respectively, L_root and L_stem are the path length from the soil to the base of the stem and the tree height, respectively, and (({phi }_{{soil}}-{phi }_{{root}})) and (({phi }_{{root}}-{phi }_{{stem}})) are the integral of conductivity from the soil to the roots and from the roots to the stem, respectively, calculated from the Kirchhoff transform.

Flow from the stem to leaves is modeled as

where a_petiole is the cross-sectional xylem area within a given petiole summed over the tree, L_petiole is the length of the petiole, (({phi }_{{stem}}-int_{0}^{L} {phi }_{{leaf}}({l}_{a})frac{d{l}_{a}}{{L}_{a}})) is the integral of the conductivity from the stem to the petiole, L_a (m² m^-2) is the leaf area index, l_a (m²) is the index of a given leaf layer, and dl_a/L_a represents the xylem area per unit leaf. Assuming there is only one leaf layer and all photosynthesis is carbon limited only, this equation simplifies to

where (({phi }_{{stem}}-{phi }_{{leaf}})) is the integral of the conductivity from the stem to the petiole under the assumption of one leaf layer. Flow from the leaf to the atmosphere is modeled as

where a_leaf is leaf area, g_s (mmol m⁻² s⁻¹) is stomatal conductance, and D (Pa) is VPD. Stomatal conductance, g_s, is modeled following ref. ⁶⁷. as

In this equation, C_a (ppm) is the atmospheric CO₂ concentration; c₁ (Pa), D₀ (Pa), and Γ (ppm) are empirical constants from the Leuning model⁷⁸; A_n (kg C month⁻¹) is net photosynthesis; and ψ_leaf is leaf water potential. The function β(ψ_leaf) serves to down-regulate photosynthesis under water stressed conditions and is determined by the carbon cost of sustaining negative water potential and loss of conductivity in the xylem. For simplicity, we assumed that β(ψ_leaf) varies linearly with the Kirchhoff transform as

where ϕ_max is the integral of maximum hydraulic conductivity of the xylem. β(ψ_leaf) varies between 1 (leaf at full hydration) and 0 (leaf under full water stress). The denominator ϕ_max is defined in terms of the maximum hydraulic conductivity (K_max) as follows:

where K_max is the model equivalent of the maximum whole-plant hydraulic conductivity (K_w,max) and b1 (% MPa⁻¹) and b2 (MPa) are the slope of the percent loss of the conductivity curve and the pressure at which 50% of xylem function is lost, respectively. Here, β broadly conforms to the solution to the Leuning model, but with a more mechanistic representation of soil moisture stress through soil water potential’s effect on leaf water potential.

The method of solution is the same as in Trugman et al.³⁹. In this way, computation of A_net is related to three climatic variables (Ψ, VPD, and CO₂ concentration), dimensions of the water conducting tissue of the tree, and tree physiological parameters.

We modified the Trugman et al. model to include a tree-liana pair and to improve the realism of the relationship between climate and plant water flow. In contrast to the use of this model for computing A_net as in Trugman et al., we use the model to define K_w,max(req), the required maximum whole-plant hydraulic conductivity, by iteratively finding the minimum K_max (Eq. 5) to yield a positive A_net on an annual timestep (Methods: Simulations). To emphasize the independence of the maximum hydraulic conductivity in the model (K_max) from plant branch-level measurements and differentiate this term in the model from K_s,max (observed branch hydraulic conductivity), we designate this term maximum whole-plant hydraulic conductivity (K_w,max) hereafter. The hydraulic conductivity variables we consider in this manuscript (K_s,max, K_w,max, and K_w,max(req)) are defined in Supplementary Table 6.

We modified the model to account for the liana growth form in three ways: inclusion of liana-tree competition, development of liana-­specific allometry, and development of a turnover routine. Our model assumes the liana and the tree are in direct competition for light and soil water. The liana-­tree pair was assigned a total leaf area of 200 m² and we varied the proportion of the total leaf area given to each the tree and the liana (Methods: Simulations). Tree and liana leaves are distributed homogeneously through the canopy and the model assumes all leaves are sun leaves. Light competition is only dependent on the quantity of leaves apportioned to each growth form; the placement of the leaves is not considered. The growth forms compete for soil water by extending a fine root area proportional to leaf area into a single, homogeneous soil layer. There is assumed to be no parasitic effect of the liana on the tree.

Liana stem length does not depend on diameter at breast height (DBH), consistent with previous modeling efforts⁵⁵. Instead, we assume liana length is as long as tree height, therefore making their canopies of equal height⁵⁵. Liana stem length may be substantially longer than tree height⁴⁷; our estimates of K_w,max(req) should be interpreted as conservative estimates. Liana DBH is then treated one of two ways. In Fig. 2, we investigate the simultaneous effects of allometry (i.e., Huber value) and hydroclimate on K_w,max(req). In this figure, we defined the total leaf area shared by the tree and the liana (200 m²) and allowed liana DBH to vary between the minimum and maximum liana DBH (1.86 and 10.7 cm, respectively) observed during a field survey in Guanacaste, Costa Rica. We then computed Huber value by dividing the sapwood area (a function of DBH) by the total leaf area apportioned to each growth form. In all other model simulations, we assigned liana DBH according to the competition scenario: 2.65 cm for the “established” scenario (equal to the mean of the observations from Guanacaste, Costa Rica) and 2.00 cm for the “invasion” scenario (the minimum stem diameter for a canopy liana; see “Model parameterization”).

We developed a turnover routine to account for differences in leaf and stem turnover between trees and lianas. The routine works as follows: during a given month, a small amount of stem is lost from an initial stem length at the beginning of the year (model parameter Lx), which corresponds with one-twelfth of the average annual stem turnover of the tree or liana. If net primary production (NPP) is negative for the month, all leaves are dropped (leaf area = 0) for the growth form and net primary production (NPP) is recalculated with leaf respiration = 0. If NPP is still negative after leaves are dropped, then stem turnover is increased to simulate a water stress response, which reduces stem respiration. This routine serves two purposes. First, the leaf turnover component allows us to account for the possibility of different phenological strategies between growth forms^79,80. To the extent possible, we allow lianas to retain leaves during the dry season to account for the potential of a “dry season growth advantage,” during which lianas maintain photosynthesis under drier conditions than trees²⁵. Second, the stem turnover component represents the fact that lianas are documented to have more rapid woody turnover than trees⁴.

The second part of our model modification is the downscaling of the model to an hourly step. The original model took as inputs VPD and Ψ at a monthly timestep. However, this does not account for the strong subdaily variation in VPD. Therefore, we modified the hydroclimate drivers of the model to account for hourly variability in VPD: Ψ remained a vector of monthly averages, while VPD became a matrix of hourly x monthly values. For use in the model, a moving average of VPD with the previous hour’s VPD was calculated to smooth the effect of increasing VPD during the day and to account for our specification of 6:00-­18:00 as daylight throughout the year.

We downscaled by computing respiration and gross primary production (GPP) for each hour of the day. GPP was set to 0 during the night (18:00-­6:00) to produce a 12-h light-­dark cycle. We summed hourly respiration and GPP to produce daily and monthly values. Then, respiration and GPP entered the turnover routine. Finally, net primary production (NPP) was computed as NPP = GPP -­ respiration.

Model parameterization

The only model inputs that differed between the tree and liana growth forms were maximum whole-­plant stem-­specific hydraulic conductivity (K_w,max mmol m⁻¹ s⁻¹ MPa⁻¹), DBH (cm), leaf area (m²), turnover (% year⁻¹), and initial stem length (m) (Supplementary Table 7). We chose to keep P₅₀ and the slope of the percent loss of conductivity (PLC) curve (model parameters b2 and b1, respectively, Supplementary Table 7) the same between growth forms because (1) our meta-­analysis suggested that the difference between growth forms in these traits is minimal compared to K_s,max (Fig. 1), and (2) this decision minimized the number of parameters contributing to differences in required K_w,max (K_w,max(req)) between growth forms.

We tested for correlations among the three traits within our meta-­analysis. We found only weak correlations between K_s,max and P₅₀ and between K_s,max and slope of the PLC curve (both: R² ≈ 0.1, p < 0.05, Supplementary Fig. 5), suggesting that fixed values for P₅₀ and slope of the PLC curve are appropriate for our analysis. Meanwhile, the correlation between P₅₀ and slope of the PLC curve is strong (R² ≈ 0.7, p < 0.05, Supplementary Fig. 5), reinforcing the fact that assigning values for these parameters with a fixed relationship best represents plant physiology. We therefore pooled observations of slope of the PLC curve and P₅₀ from both growth forms in our meta-­analysis to compute our estimates of b1 and b2.

DBH distributions and average DBH for each growth form were taken from surveys of lianas and trees in a second growth forest of Guanacaste, Costa Rica (Supplementary Fig. 8). For the scenario of an established liana in a tree canopy (“established scenario”), we assumed a liana DBH equal to the mean observed at Guanacaste, ≈2.65 cm. For the scenario of a liana invading a tree canopy (“invasion scenario” considered in the Supplementary Discussion), we assumed a liana DBH of 2 cm⁸¹. In all simulations, tree DBH was assumed to be the average from the tree survey of Guanacaste (≈18 cm). For the established scenario, we assumed the liana occupied 40% of the total leaf area (80 m²) and the tree occupied 60% of the total leaf area (120 m²). For the invasion scenario, we assumed the liana occupied 10% of the 200 m² total leaf area (20 m²) and the tree occupied 90% of the total leaf area (180 m²).

For most traits, there was limited evidence for tree-­liana differences (e.g., wood density, Glass’ ∆ < 1) or there were insufficient data to parameterize the liana growth form (e.g., root:shoot ratio). Specific leaf area (SLA) was a special case. Although SLA was found to be significantly different between growth forms (Glass’ ∆ ≈ 1), we did not assign different values to lianas and trees because the TRY results are likely influenced by the low SLA of desert-­dwelling and montane shrubs within the tree growth form. Values of the inputs and parameters that differ from the original model are provided in Supplementary Table 7^{3,82,83,84,85,86}. All other parameters are the same as those used in the original model³⁹.

Sensitivity analysis

We conducted an extensive sensitivity analysis of our model to identify the parameters that are most influential to determining K_w,max(req). For each parameter in the model (n = 25), we computed K_w,max(req) with a 50% reduction and 50% increase from the default (mean) value while holding all other parameters at their default values. We then found the difference between K_w,max(req) from the 50% increase and 50% decrease in parameter value and divided the difference by the K_w,max(req) at the default parameter value; we report this computation as the “sensitivity.” We computed the sensitivity of each parameter for two hydroclimate conditions, BCI and Horizontes, and for the two competition scenarios, established and invasion (with respect to the liana). When tree K_w,max(req) was computed, we held all liana parameters at their default values and vice versa. This amounted to over 400,000 additional annual model simulations. This sensitivity analysis informed the parameters that we used field collected data to constrain, including diameter at breast height (DBH), P₅₀, and hydraulic path length (Lx). Where constraining the parameters with field data was not possible, we conducted additional model simulations with alternative scenarios. For example, given that we found the model to be sensitive to the total leaf area, we ran additional simulations to create Fig. 3 and 4 under alternative total leaf areas, 150 m² and 400 m². The Supplementary Methods (Supplementary Method: Sensitivity analysis) offers more detailed results of our sensitivity analysis and how those results informed our modeling procedure.

Climate data

Our model computes NPP as a function of carbon dioxide concentration ([CO₂]), Ψ, and VPD. We set [CO₂] at 400 ppm, a low-end estimate for the 21^st century, for all model simulations except in our predictions of 2100, in which [CO₂] = 550 ppm⁸⁷. Our hydroclimate data come from two Neotropical forests with contrasting hydroclimate conditions, Horizontes, Costa Rica and Barro Colorado Island (BCI), Panama. Ψ was determined from Medvigy et al.⁸⁸. (Horizontes) and Levy-­Varon et al.⁸⁹. (BCI). In each case, Ψ was estimated for multiple soil layers in the original dataset. However, because measurements were not taken at the same soil depths at each location and because our model assumes there is only one soil layer, we used Ψ estimates from only the 15 cm depth, which was available for both sites, for all simulations. VPD data are from a reanalysis data product for Horizontes averaged over 2007-­2018 (ref. ⁹⁰). For BCI, data are from the Lutz Tower from 27 May 2002 to 5 June 2020 at 48 m canopy height. We computed VPD from relative humidity and air temperature data at both sites as follows:

where SVP is saturation vapor pressure (hPa), AT is air temperature (°C), RH is relative humidity (%), and VPD is in Pa. At both sites, VPD data were averaged across year and day of the month. Changes in monthly Ψ and VPD for BCI and Horizontes are available in Supplementary Fig. 9.

Simulations

We used our model to simulate required conductivity (K_w,max(req)) by identifying the smallest value of whole-­plant conductivity (K_w,max) at which NPP is positive under the given hydroclimate and liana-­tree competition conditions. To compute K_w,max(req) we performed the following steps: (1) define the simulation inputs, including DBH and total leaf area fraction for each growth form, and hydroclimate (i.e., VPD and Ψ); (2) run the model for each month with the smallest value of K_w,max available; (3) sum the monthly NPP computed by the model; (4) if NPP > 0, define K_w,max(req) as the current K_w,max; and (5) if NPP ≤ 0, select the next lowest value of K_w,max and repeat the steps until NPP > 0, at which point K_w,max(req) is identified (Supplementary Fig. 2).

We emphasize that the model depends on K_w,max, whereas it is much more common to measure terminal branch K_s,max. Because of the uncertainty associated with scaling between K_w,max and K_s,max, our estimates of K_w,max(req) should be compared to observed K_s,max with caution. To reduce uncertainty in this parameter, we urge further measurements of K_w,max.

We first simulated K_w,max(req) under different hydroclimate scenarios, as shown in Fig. 2. The hydroclimate scenarios are tropical dry forest and tropical moist forest (Methods: Climate data). Instead of defining the liana DBH, we computed K_w,max(req) over a range of DBH values observed in our liana survey dataset from Horizontes, which allowed us to avoid assigning a fixed allometry to lianas in our initial simulations.

We similarly computed K_w,max(req) under a variety of VPD-­Ψ scenarios. The indices were computed by linearly interpolating the hydroclimate between the driest (Horizontes) and the wettest (BCI) sites for a length of 100 indices for both VPD and Ψ. For each combination of Ψ and VPD (10,000 combinations), we computed K_w,max(req) using the method outlined above.

We extended our computation of K_w,max(req) for each growth form into the future under a gradient of increasing VPD conditions. Because of the high uncertainty surrounding the magnitude of increases in VPD over the next 100 years¹², we computed K_w,max(req) under a variety of VPD scenarios, ranging from 10% to 100% increase in VPD from the present at Horizontes (Supplementary Fig. 3). For the model simulations involving future VPD scenarios, we additionally changed the atmospheric [CO₂] to 550 ppm to reflect the dependence of climate change (i.e., increasing VPD) on increasing atmospheric [CO₂].

Finally, we investigated the potential influence of liana and tree physiological adaptation to drying hydroclimate via adapting P₅₀. Because of the strong empirical correlation between P₅₀ and the slope of the percent loss of conductivity curve (Slope), we simultaneously varied these two parameters in three additional scenarios, with hydroclimate conditions predicted for 2100 (i.e., 100% increase in VPD, no change in Ψ, [CO₂] = 550 ppm). We used the “established” competition scenario and assumed the same adaptation scenarios for both liana and tree K_w,max(req) simulations. The three scenarios are as follows: b1 = 0.92% MPa⁻¹, b2 = −2.25 MPa; b1 = 0.73% MPa⁻¹, b2 = −2.5 MPa; and b1 = 0.49% MPa⁻¹, b2 = −3 MPa. The most extreme liana P₅₀ observed in the literature we included in our extended meta-analysis is -2.99 MPa; therefore, our P₅₀ adaptation scenarios are consistent with the most drought-tolerant observations of present-day liana P₅₀.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com