Basin-wide variation in tree hydraulic safety margins predicts the carbon balance of Amazon forests – Nature

Site description

We assemble a pan-Amazon dataset of key HTs ((varPsi )₅₀, HSM₅₀ and (varPsi )_dry), including 129 species distributed across 11 forest sites (Fig. 1 and Extended Data Fig. 1). The sites are old-growth lowland forests (less than 400 m of elevation), with no evidence of significant human disturbance, located in western, central eastern and southern Amazon. They were specifically chosen to span the full Amazonian precipitation gradient and to encompass the principal axes of species composition in the Amazon. The MAP varied from around 1,390 to around 3,170 mm yr⁻¹ and mean MCWD varied from −640 to −15 mm across sites. Summary information for all sites can be found in Supplementary Tables 1 and 2.

Species selection

To characterize drought sensitivity across a wide set of species and strategies, we sampled the most dominant adult canopy and subcanopy tree species at each site. For TAP, MAN and CAX, we used published data from refs. ^6,27,30 which follow the same methodology as this study. The sampling effort at each site varied from 7 to 26 species which represented between 14% and 70% of the total basal area (Supplementary Table 2). Sites for which less than 30% of the total basal area was sampled (ALP-1, ALP-2, SUC, CAX and MAN) are hyperdiverse forests and lack the clear dominance structure by a few species observed in less diverse plots (for example, in the southern Amazon NVX site, the seven species sampled account for more than 50% of the basal area). Previous work by ref. ⁶, show that the MAN site, despite having the lowest sampled basal area of all sites presented in this study (about 14%) is representative of the broader floristic community, as adding a broader array of species-level hydraulic trait data did not significantly change basal area weighted mean (CWM) values. The same study found that mean species values are not likely to differ from community mean values if (1) species dominance is not driven by a few species, (2) traits have low dispersion around the mean (low standard deviation compared to the mean) and (3) traits are randomly distributed across species dominance distributions. For the other four sites for which sampled coverage was less than 30%, these criteria are generally satisfied (for example, cumulative dominance of the five most dominant species at ALP-1 is 27.9%, ALP-1 26.2%, SUC 15.0% and CAX 10.7%, standard deviation of (varPsi )₅₀ is between 32% and 49% of the mean value at each site and there is no relationship between species dominance and HT. Thus, basal area weighted mean trait values for the 11 sites probably well represent the broader unsampled community of trees.

Abiotic data

To characterize climatological water deficit at each site, we calculated the MCWD⁶⁹, which is a widely used measure of dry season intensity for Amazon forests^13,16,70 that expresses the cumulative water stress experienced within an average year⁶⁹. The MCWD metric assumes that a forest experiences water deficit if monthly precipitation does not meet evapotranspirational requirements and accumulates that deficit over all successive months with rainfall lower than evapotranspiration (E) values⁶⁹. Monthly water deficit (WD_n) was then calculated as the difference between precipitation (P) and evapotranspiration demand in each month n. MCWD was computed as the maximum monthly cumulative water deficit (CWD) experienced over an average year, for which the change in water deficit in any given month n is calculated as the difference between precipitation falling that month (P_n) and an assumed evapotranspiration demand (E_n, mm month⁻¹). For any given month n,

As all of our plots are in the southern hemisphere, their hydrological year coincides with the calendar year, allowing us to start our MCWD calculations at the beginning of each calendar year. For statistical analyses, we use the long-term mean MCWD for each location. Monthly precipitation data were obtained from the tropical rainfall measuring mission (TRMM TMPA/3B43 v.7)⁶⁴ at 0.25° spatial resolution from 1998 to 2016. To estimate evapotranspiration, we used monthly ERA-5-Land Reanalysis E data at 0.1° spatial resolution from 1998 to 2016⁷¹, as this product has been suggested to well represent evapotranspiration estimates in the Amazon⁷². To have one value of evapotranspiration demand per site (E_n in equation (1)), we used the mean E value for the 3 months with highest E across years. Mean annual temperature data at 1 km spatial resolution were obtained from Worldclim2 (ref. ⁷³).

We performed an alternative assessment computing MCWD on the basis of MOD16 (ref. ⁷⁴) evapotranspiration product and on E estimation of 100 mm per month⁶⁹ and we also computed MAP on the basis of TRMM⁶⁴ and CRU⁷⁵ data. The main results remained similar, independent of the climate product used (Supplementary Table 7).

Collection of plant material

One fully sun-exposed top-canopy branch (or branch at the maximum height reachable by climbers) was collected from, on average, three individuals of each species at each site for subsequent construction of xylem vulnerability curves. The same or a second set of branches, in the same canopy position, was used to extract samples of wood density and LMA. For embolism resistance determination, data collection was undertaken during the wet season, when forests were maximally hydrated. Branches (more than 1 m long) were harvested during predawn or very early in the morning, to capture a fully hydrated starting point. Immediately after collection, basal portions of branches were wrapped with a wet cloth and branches were placed in a humidified opaque plastic bag to avoid desiccation during transport. Bags were sealed and carried to the field station for determination of xylem vulnerability curves. For samples not collected during predawn, branches were placed in a bucket, recut under water, covered with an opaque plastic bag and left to rehydrate for at least 5 h.

Xylem embolism resistance ((varPsi )
₅₀ and (varPsi )
₈₈)

To quantify xylem resistance to embolism of Amazonian trees species, we focused on the water potentials associated with (varPsi )₅₀, given its wide use as a critical embolism resistance threshold^4,5. To derive this parameter, we constructed xylem vulnerability curves by simultaneously measuring percentage of embolism formation and xylem water potential under progressive desiccation⁷⁶. We estimated embolism using the pneumatic method of ref. ⁷⁷, which quantifies the air extracted from within branches at each stage of dehydration and expresses this as a percentage of air discharge (PAD), defined as the percentage difference between the maximum amount of air removed under extreme dehydration (100% PAD) and the minimum amount removed under maximum hydration (0% PAD)⁷⁸. For our measurements we used manual, self-constructed pneumatic devices, following ref. ⁷⁸. Although automated devices for measuring air discharge are now available, these were not available at the time of our data collection. For all air discharge determinations, we applied the protocol of ref. ⁷⁷ whereby measurements of air discharge were made over a 2.5 min interval. We note that the absolute volumes of air discharged are sensitive to the time interval of the discharge measurements, as shown by ref. ⁷⁹,who report a difference of about 10% on the absolute air discharge measured for 15 s versus 115 s. There are still methodological uncertainties that require further investigation, including how the contribution of extraxylary discharge varies across different Amazonian species. Recent work using a pipe pneumatic model to simulate gas diffusion from intact conduits suggests that the overriding source of discharged air is from embolized xylem vessels although there is a small contribution (estimated to be about 9% over 15 s of discharge) from extraxylary pathways⁸⁰. It is also important to note that the method measures embolism from vessels connected to the cut end of the branch from which gas is sampled and that there may be more embolism from vessels that are not directly connected to the cut end⁸⁰. However, embolism spread during the branch dehydration method for embolism induction used in this study is expected to be predominantly from the cut branches⁷⁶ and is corroborated by the strong agreement between petiole embolism status using the pneumatic method and leaf vein embolism assessed using optical approaches⁸¹.

The portability, ease of use and low cost of the pneumatic method make it ideally suited for use in remote tropical environments in which laboratory infrastructure is often minimal. Several studies have shown that (varPsi )₅₀ values derived from the pneumatic approach agree closely with those derived using more laborious methods^{77,79,81,82,83,84}. For the TAP site in this study, (varPsi )₅₀ determinations based on the pneumatic method were compared with values derived from xylem vulnerability curves of percentage loss of conductance (PLC) constructed using a hydraulic ultralow flow meter⁸⁵ and found a strong agreement (R² = 0.83) between both methods⁸⁴, further corroborating findings from previous studies (refs. ^27,84 provide detailed description of the hydraulic method used). Although one study⁸⁶ (but see refs. ^84,87) proposed that the method may be unsuitable for long-vesseled species, we find no evidence of any vessel length bias in our (varPsi )₅₀ estimates derived from the pneumatic method (standard major axis (SMA) regression (varPsi )₅₀ versus maximum vessel length: P = 0.15, R² = 0.02).

The initial PAD measurement for each branch was made immediately after removing the branch from a sealed opaque plastic bag to ensure that vulnerability curves started from a maximally hydrated state. Subsequent measurements were then conducted successively throughout the dehydration process, with approximately eight to ten measurements per individual used to construct each curve. Branches were progressively dried through the bench dehydration technique⁷⁶. Between each dehydration state, branches were bagged for a minimum of 1 h to equilibrate leaf and xylem water potentials. Leaf water potential (used as a proxy for xylem water potential following equilibration) was measured with a pressure chamber (PMS 1505D and PMS 1000, PMS instruments).

We used the exponential sigmoidal function of ref. ⁸⁸ to calculate (varPsi )₅₀ for each species at each site:

where S is the slope of the curve, Ψ_x is xylem water potential (MPa) and Ψ₅₀ is Ψ_x corresponding to a PAD of 50%.

Following ref. ⁸⁹, we computed Ψ₈₈ as:

({boldsymbol{Psi }})
_dry and HSMs

To calculate how close Amazonian trees operate to critical embolism thresholds in nature, we measured in situ midday leaf water potentials during the peak of the dry season ((varPsi )_dry). Sampling campaigns closely corresponded with the time of most intense water deficit (Extended Data Fig. 2) and the year of sampling was not climatologically anomalous. We sampled three to six top-canopy fully expanded and sun-exposed leaves per individual (three individuals per species for 129 species in total across 11 sites) from 11:00 to 14:30. Parameter (varPsi )_dry was measured with a pressure chamber (PMS 1505D and PMS 1000, PMS instruments) in situ immediately postsampling and the values of different leaves averaged per individual. In our protocol we tried to minimize the time spent between branch cutting and the leaf water potential measurement with the pressure chamber (around 3–5 min). We collected branches (40–60 cm in length, depending on the species and leaf size) that were fully exposed to light from the top part of the canopy (highest part that the climbers could reach), from apparently healthy and undamaged individuals. Telescopic shears (normally four to six poles, with total length of 5–7 m) were used to access and cut the branches. As soon as the branches hit the ground, the branches were bagged in a black and opaque plastic bag and transported to the pressure chamber, which was located inside the plot. We then collected three to six healthy and fully expanded leaves for each individual and immediately (after the cut) placed them into the pressure chamber. All of the processes were made as quickly as possible to avoid dehydration.

Because of pressure drops in transpiring leaves, we note that the water potentials measured are probably lower than the branch water potential values at the time of measurement. Apart from aseasonal ever-wet forests, which have no climatological dry season (monthly_precip < 100 mm), data collection took place in the peak of dry season (Extended Data Fig. 2) during what were climatically normal years. For each species at each site we calculated the HSM with respect to ({varPsi }_{50}) (HSM₅₀), as the difference between species-level (varPsi )_dry, taken as the minimum (varPsi )_dry value of all individuals for that species and (varPsi )₅₀. All (varPsi )_dry measurements were made in climatologically normal years. (See Supplementary Table 8 for further information of sampling dates for each site). We also calculated the HSM with respect to ({varPsi }_{88}) (HSM₈₈) for all sampled species.

Wood density and leaf mass per area

We combined published and new field measurements of LMA and wood density to understand the power of these traits relative to HTs in predicting forest carbon balance. Stem wood density data (WD_stem) were obtained from the Global Wood Density database^65,66 and calculated as species mean values. We measured wood density at branch level (WD_branch) using a water displacement method⁹⁰. In this method, branch segments of about 25 mm length and 12 mm diameter were first cut and debarked. Samples were then placed in a recipient with filtered water to rehydrate for 24 h and subsequently weighed with a three-decimal scale. After this, the sample was oven-dried for 48–72 h at 70 °C and the dry weight measured with a balance. Wood density was then expressed as the ratio of wood dry mass and wood fresh volume (g cm⁻³). Branch wood density measurements were made in all sites except NVX. For this site, we used stem wood density values^65,66 for each of our target species.

We measured LMA for all sampled species in each of the 11 sites. For this, all leaves were detached from a selected branch and a subsample of 10–20 leaves per branch were taken, numbered and scanned. All the other leaves were kept separate to be oven-dried. This was usually done as soon as possible after returning to the field station. When it was not possible to scan the leaves straight away, we placed all the detached leaves into a sealed plastic bag in the dark and stored them for no more than 24 h. After scanning, all leaves were oven-dried for 48–72 h at around 70 °C. Once dry, the subsampled numbered leaves were individually weighed and the non-numbered leaves were weighed together with a precision scale (three decimals). On the basis of the relationship between the fresh area and dry weight of individual leaves (from the subsampled 10–20 leaves) and having the dry weight of all the leaves of the branch, we estimated the fresh leaf area corresponding to the entire branch. The LMA was then calculated as the ratio of leaf dry mass to fresh area, expressed in g m⁻². We then calculated basal area weighted mean values for all these traits for each site (Supplementary Table 2). The number of species sampled for each trait is shown in the Supplementary Table 9. Further leaf habit information of sampled species is provided in Supplementary Table 10.

Water deficit affiliation

To describe Amazonian species-level biogeographical distributions, we used published WDA data⁴⁵, which describes the spatial association of Amazonian tree species with climatological water availability. WDA was calculated as the mean climatological water deficit across inventory plots in which a species occurs weighted by its relative abundance in each of 513 forest plots broadly distributed in the western Neotropics⁴⁵. More negative WDA values represent dry-affiliated species, whereas wet-affiliated species are represented by less negative WDA values.

Forest dynamics data

We used long-term forest plots from the RAINFOR network³⁹ to help understand the relationship between hydraulic attributes and stand-scale carbon dynamics. Thus, we computed AGB net change (ΔAGB, Mg ha⁻¹ yr⁻¹), annual aboveground wood production (AGWP, Mg ha⁻¹ yr⁻¹), annual AGB mortality (AGB_MORT, Mg ha⁻¹ yr⁻¹), annual instantaneous stem mortality rate (% y⁻¹) and woody biomass residence time (τ_w) for the same forest plots sampled directly for HTs. For the two plots (MAN and TAP), for which we did not have access to forest dynamics data, we used information from a permanent RAINFOR network forest plot in the same landscape, with the most similar structure and species composition (BNT-01 and TAP-02, respectively; Supplementary Table 6) to our sampling plots. For CAX, we used published data by ref. ⁹¹ for the control plot. The other six plots are part of the RAINFOR network³⁹, having been established by and/or monitored by RAINFOR partners (Supplementary Table 5). Plot data for these analyses were curated and obtained via the ForestPlots.net database^92,93, for which standard quality control procedures are applied. We only included plots in the analysis that lacked a history of recent anthropogenic disturbance. For all forest dynamics analyses we excluded KEN plots because of a fire event that occurred in the region in 2004⁶⁸ and may still be affecting biomass stocks and dynamics. Following previous studies^15,94, plots smaller than 0.5 ha that were up to 1 km apart from each other were combined and treated as a single plot (for example, TAP-54, TAP-55, TAP-56 and TAP-57 treated as TAP-02, the plot we used to represent TAP). For each plot, we only included pre-2015 El Nino censuses and selected the census start date to be as consistent as possible across plots. For this we excluded pre-2000 measurements, apart from TAP plot for which censuses were available only from 1983 to 1995. For other plots, we used the earliest census available for this plot if data collection started after 2000 (VCR-02 plot, for example, which starts in 2003). We tried to ensure that biomass dynamics metrics used in the analyses represented at least 10 yr of total monitoring time per plot. If application of the 2000 start date for a given plot resulted in fewer than 10 yr of monitoring, we also included the census date immediately before 2000 (99 for BNT-02 plot, which we used to represent MAN) to ensure at least 10 yr of monitoring (Supplementary Table 5). The monitoring time used for the plots included in the analysis was on average 12.3 (s.d. = 2.5) yr. In RAINFOR plots, all live individuals of more than 10 cm in diameter at breast height (DBH) are repeatedly measured over time, using standardized protocols, with species identified and careful records kept of trees that die or recruit from one census to the next. AGB for each census per plot was computed using the ref. ⁶³ equation for moist forests on the basis of tree diameter, wood density and height. As local height data were often unavailable, a Weibull equation with regionally varying coefficients was used to estimate height following ref. ¹¹. Species-level wood density values from the Global Wood Density database^65,66 were used to compute AGB, AGWP and AGB_MORT. For each census, biomass values were calculated for all dicotyledonous trees in the plots above the 10 cm DBH cut-off and summed to give total stand-level biomass stocks.

We estimated annual ΔAGB (Mg ha⁻¹ yr⁻¹) for a given plot as the difference in AGB between the final and initial census used (AGB_{final census} − AGB_{initial census}) divided by the monitoring length (Date_{final census} − Date_{initial census}) in years. For each census interval per plot, we also computed annual AGWP (Mg ha⁻¹ yr⁻¹), following ref. ⁹⁵, which encompasses (1) the sum of the growth of surviving trees, (2) the sum of AGB of new recruits, (3) the estimated sum of growth of unobserved recruits that dies and (4) the estimated sum of unobserved growth of initial trees that died, within a plot in a given census interval, divided by the census interval length (yr) (see also ref. ⁹⁴). For each plot, we computed annual AGB_MORT, including unobserved components, which is defined as the sum of the AGB of all dead trees, plus the estimated growth of recruits that died before they could be recorded in the second census and the sum of estimated unobserved growth of trees that died within an interval, divided by the census interval length⁹⁴.

As AGB varies across sites it is useful to account for this when comparing sites. We therefore also computed relative ΔAGB (ΔAGB/AGB), relative AGWP and relative AGB_MORT by dividing absolute values by the time-weighted mean standing woody biomass across censuses per plot. Both absolute and relative values are presented in the Extended Data Figs. 7, 8 and 10 and Supplementary Table 4). We computed the annual instantaneous stem mortality rate (% y⁻¹) following ref. ⁶⁷:

in which A is the number of stems per ha in the beginning of the census interval and B is the number of stems per ha that survived throughout the census interval. Owing to the sensitivity of these rates to census interval effects, we standardized them to a common census interval, following ref. ⁹⁶. For all calculations above (AGB, AGWP, AGB_MORT and stem mortality) we used the BiomasaFP R package⁹⁷. We calculated the time-weighted mean values of all these absolute and relative parameters (AGB, AGWP, AGB_MORT and stem mortality) to have one value per plot. We then calculated woody biomass residence time (τ_w) as the ratio of the time-weighted mean standing woody biomass and the time-weighted mean annual biomass mortality⁵².

To test whether relationships between HSM₅₀ and forest dynamics at plot level apply over landscape scales and to account for the influence of within- and among-plot stochasticity in dynamics, we also we used mean forest values of forest dynamics metrics across groups of plots (clusters) in the same landscape with similar structure and composition to plots sampled for hydraulic measurements (Supplementary Tables 5 and 6). For this cluster-level analysis, we excluded white-sand forests and permanently water-logged swamp forests because they are extreme edaphic habitats, known to have a more limited and edaphically specialized tree flora⁹⁸. We also excluded forests lying within active floodplains of rivers because their flora is also distinctive and, like swamp forests, they have access to more water beyond that which is climatically determined. In total, we used data from 34 long-term monitoring plots (31.37 ha of forest). For this analysis, we used cluster mean forest dynamic values (instead of plot cluster weighted mean, for example) because plot area and monitoring length did not vary considerably within clusters (Supplementary Table 5). To account for sampling effort variation across cluster of forest plot, we tested if the residuals of the relationship between relative ΔAGB and HSM₅₀ were related to cluster mean monitoring time (mean ± s.d. was 12.1 ± 1.8 yr) and cluster total area (3.9 ± 3.0 ha). No weights were assigned to each data point in the regression because we found no evidence of relationships between the residuals and sampling effort across clusters.

Statistical analysis

To examine the distribution of HTs ((varPsi )₅₀, (varPsi )_dry and HSM₅₀) across Amazonian tree taxa (N = 129 species), trait values were averaged for species occurring at several sites. We conducted statistical analyses to investigate differences in species-level hydraulic trait values among different forest types and geographical regions and also to evaluate controls of water availability on basal area weighted mean HT across the study sites. To examine differences in HTs among forest types, we first grouped our 11 forest sites into three forest types, based on DSL: (1) ecotonal long DSL forests—DSL equal to 6 months, MAP and MCWD less than 1,600 and −470 mm, respectively; (2) intermediate DSL forests—DSL ranging from 5 to 2 months, MAP between 1,990 and 2,650 mm and MCWD varying from −288 to −184 mm; and (3) ever-wet aseasonal forests—DSL about 0 months, MAP and MCWD greater than 2,950 and −15 mm, respectively (Extended Data Fig. 1 and Supplementary Table 1). To test for statistical differences in HTs across forest types, we performed a one-way Kruskal–Wallis followed by a post hoc Mann–Whitney–Wilcoxon rank sum test. Western and central eastern Amazon forests have fundamentally different dynamics in that western Amazon forests are characterized by high growth and turnover whereas central eastern forests are associated with slow growth and turnover^35,36. To test for differences between species in intermediate DSL sites in western Amazon (FEC and TAM) and central eastern Amazon (CAX, MAN and TAP), we performed Wilcoxon rank sum tests. Linear models were constructed to evaluate relationships between basal area weighted mean HT and MCWD (Supplementary Table 7). For all analyses, we use a significance level of 0.05.

To investigate if species biogeographical distributions are related to mean HT, we used SMA regressions with WDA as the response variable. Following Esquivel-Muelbert et al.⁴⁵, we restricted our analysis to the western Amazon as these published WDA data are based entirely on species distributions within western Amazon, helping to control for the potentially confounding effects of differences in soil and forest dynamics across Amazonian regions. Our subsample for this analysis encompassed a total of 87 species distributed across aseasonal, intermediate DSL and ecotonal long DSL forests, with MAP across plots ranging from 1,390 to 3,170 mm. SMA regressions were performed using the smatr package⁹⁹ in R.

Using our entire dataset across the Amazon, we evaluated whether HTs were better predictors of Amazon forest carbon balance than climatic factors or other leaf and wood traits. More specifically, we performed bivariate SMA models to investigate relationships between HTs ((varPsi )₅₀, (varPsi )_dry and HSM₅₀), climatic data (MCWD, MAP, DSL and MAT) and other functional traits (LMA, WD_stem and WD_branch) versus long-term ΔAGB at plot level. We computed basal area weighted mean LMA, WD_branch and WD_stem data^65,66. To account for the influence of multiple testing, we applied a Bonferroni correction to P values for bivariate regressions. SMA models were further conducted to examine the relationship between HSM₅₀ versus absolute and relative values of AGB annual woody production, AGB annual mortality, stem mortality and residence time of woody biomass. Supplementary Table 5 presents summary information per plot and clusters. All presented analyses were performed in RStudio v.1.1.423 (ref. ¹⁰⁰).

Reporting summary

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