Rainfall and topographic position determine tree embolism resistance in Amazônia and Cerrado sites

Droughts are predicted to increase in both frequency and intensity by the end of the 21st century, but ecosystem response is not expected to be uniform across landscapes. Here we assess the importance of the hill-to-valley hydrologic gradient in shaping vegetation embolism resistance under different rainfall regimes using hydraulic functional traits. We demonstrate that rainfall and hydrology modulate together the embolism resistance of tree species in different sites and topographic positions. Although buffered by stable access to groundwater, valley plants are intrinsically more vulnerable to drought-induced embolism than those on hills. In all study sites, the variability in resistance to embolism is higher on hills than on valleys, suggesting that the diversity of strategies to cope with drought is more important for tree communities on hills. When comparing our results with previously published data across the tropics, we show greater variability at the local scale than previously reported. Our results reinforce the urgent need to extend sampling efforts across rainfall regimes and topographic positions to improve the characterization of ecosystem resistance to drought at finer spatial scales.


Introduction
Understanding the vulnerability of trees to drought is critical to understand the response of tropical ecosystems to projected drier and warmer conditions [1][2][3][4]. While reduced rainfall and increased vapor pressure deficit are key drivers of tree mortality [5][6][7], topography creates spatial heterogeneity in plant groundwater access and mediates contrasting vegetation response to droughts [8][9][10][11][12]. Field observations show that trees in valleys are often less affected by droughts than their upland counterparts, with  [36][37][38]. All water entering the system is a balance between precipitation (determined by regional climate) and evapotranspiration (determined by climate and vegetation). After rainfall, water that is not intercepted by the canopy infiltrates the soil, and part of it is stored in the unsaturated zone, which encompasses the soil above the water table (WT), where soil pores are filled with air and water. The other part is stored below the WT, in the saturated zone where all the soil pores are filled with water. Water flows downhill in the subsurface through groundwater lateral flow, below the WT, keeping a shallow WT in the valleys, accessible to plants. The rise and fall of the WT tracks the rainfall seasonality with a lag [8,28,39]. Under the hills, even at its peak, it is still too deep for most plants. In the valleys, the WT is shallow year-round, and can become too shallow in the wet season, creating waterlogging. This different exposure to hydrologic stress (drought under the hills, waterlogging in the valleys) is hypothesized to select for plant species according to their adaptations to hydrologic stress [34]. evidence from the tropics [13][14][15][16][17] and temperate areas [9,18,19]. Furthermore, remote sensing-based and modeling studies have also showed that valley communities are less sensitive to drought and show lower drought mortality [20][21][22]. Contrasting evidence of higher drought mortality in valleys [23,24] and water table (WT) variability negatively affecting trees in lowland positions [25] raise questions as to whether the hill-to-valley gradient in water access is equally important in buffering plants from droughts in regions with different rainfall regimes, topography and species composition. This highlights the need to look at how rainfall, topography and species composition interact together to shape community-level response to droughts.
Hydrologic conditions at the landscape scale can be remarkably different from what is predicted by rainfall [26]. The slow motion of water in the subsurface from high to low grounds creates patches of wetter and drier terrain under the same climate (figure 1). Valleys receive water from neighboring hills mostly through surface and subsurface flow, and thus have a shallower WT even during the dry season and extreme drought events [8,15,27,28]. On the hills, the WT is deep and out of reach for most trees, which rely either on adaptations to withstand drought or mechanisms to avoid it, such as deep rooting, accessing moisture stored in the deep soil [29][30][31][32][33]. The resulting gradient in water access from hills to valleys creates hydrologic niches, which species occupy according to their adaptations to drought and waterlogging [34,35]. At the regional scale, mean annual precipitation (MAP) and rainfall seasonality influence the distribution of floristic composition and physiological traits related to drought tolerance [40][41][42][43][44][45], with lower rainfall and higher seasonality being associated with higher drought resistance [43][44][45][46]. At the local scale, hydrologic gradients drives floristic [42,47,48] and functional turnover, with hill species having more embolism-resistant xylem than valley species [36][37][38][49][50][51][52] and thus able to survive drier conditions. However, the role of the hydrologic gradient in modulating vegetation drought tolerance across different rainfall regimes remains largely understudied.
Here, we evaluate the role of hydrology in shaping embolism resistance of tree species in valley and hill positions across a gradient of decreasing rainfall and increasing seasonality. We use the P 50 metric, or the plant water potential at which 50% of water transport capacity is lost due to the formation of air bubbles in the vessels (i.e. embolism). Lower (more negative) values of P 50 have been associated with higher resistance to drought and less mortality [53][54][55]. We measured P 50 in dominant trees in valleys and hills for three sites in South America, spanning a rainfall gradient from 1370 to 2245 mm yr −1 associated with an increase in dry season length and intensity. We propose two hypotheses: (H1) Valley species have lower embolism resistance than hill species even under low precipitation and high seasonality, as the WT remains shallow and thus accessible to valley species despite a more stressful rainfall regime. (H2) The hydrologic gradient's influence on P 50 depends on local rainfall and topography. Rainfall is the first-order control on how much water is available for vegetation, while topography reorganizes species according to their drought tolerance. Thus, the dissimilarity in embolism resistance between hill and valley species increases with rainfall seasonality, as hills are progressively more stressed by reduced rainfall and soil moisture, while valleys are still buffered by a shallow WT.

Study areas and environmental variables
We tested our hypotheses at three sites distributed along a MAP and seasonality gradient (figure 2, red dots). To characterize the rainfall regime in each site, we use 30 year rainfall records from the Agência Nacional dasÁguas in each site (table S1 To monitor the position and dynamics of the WT in the Tapajós and PNCV sites, we installed three groundwater wells at the lower portion of the hillslopes using a hand auger at the end of the dry season in October 2020 in the PNCV site and November 2021 in the Tapajós site. Wells were drilled below the position of the WT, to ensure continued measurements, and were spaced by ∼20 m starting close to the stream and continuing along the hillslope. The depth to the WT was measured biweekly to monthly in the PNCV from December 2020 to November 2022 and monthly in Tapajós from November 2021 to June 2023 with a common measuring tape. Because our depth of investigation was limited by the reach of the hand auger (5-8 m), we were unable to install wells in the upper portions of the slope and the hills. To obtain WT depth under the hills, we fitted an analytic hydrological model to the well data in Tapajós and PNCV (SM1.1, table S2). For the Manaus site, we obtained hill and valley WT depth from the literature [56].

Species selection
For Manaus, we used data from [36], which sampled 28 tree species (table S3) by their topographic affiliation (hill or valley). Of those, half are among the 50 most abundant genera in Reserva Ducke and 12 are among the 227 hyperdominant species in Amazônia [40]. We also included 17 tree species sampled by [46] on hills in Reserva Cuieiras to supplement the data, given the enormous species diversity of Amazonian forests. Reserva Cuieiras is located 60 km from Manaus and is in the same rainfall envelope, also sharing similar geology and soils [46]. We therefore group species from Reserva Cuieiras and Reserva Ducke as representative of one site (Manaus). Species were selected according to their dominance [46] and/or topographic affiliation [36].
In Tapajós, we used sampling plots from one of the permanent research modules previously installed by the Biodiversity Research Program (PPBio Santarém, https://ppbio.inpa.gov.br/nregionais/ nrsantarem). Plots measured 1 ha and followed the PPBio/RAPELD methodology [57], which establishes plots along contour lines (thus controlling topography) and is regarded as a standard sampling methodology for long-term research sites in Amazônia, such as those sampled in Reserva Ducke. Ten plots were installed in the km 117 A module, situated 50 km south of the long-studied Large-Scale Biosphere-Atmosphere Experiment in Amazônia (LBA) research site. We selected species from two plots (one hill, one valley) based on relative dominance in each community, which gave us 13 tree species in each plot (table S3), corresponding to 31.9% and 29.3% of the community basal area in the valley and hill, respectively.
In the PNCV site, we established plots on two positions along a hillslope. The hill plot is in a welldrained Cerrado sensu stricto formation and the valley plot is in a narrow riparian forest bordering a perennial stream. Because of the short length of the hillslope and the narrow nature of Cerrado riparian forests [58], we established 4 plots of 5 m × 25 m, separated by 25 m, in each location, with the long axis following the topographic contour. We surveyed every individual tree with a diameter at base ⩾2 cm (hill) or diameter at breast height ⩾10 cm (valley). We chose the five most dominant tree species in both the hill and valley plots (table S3), representing 75% and 67.4% of the community basal area, respectively.
While we adopted the same procedure for species selection on our sites, differences in species richness and dominance lead to different basal area representation. The high species richness and homogeneous dominance in Amazonian forests [59] would require a large number of species to obtain the same representation as in the PNCV, where tree species richness is low and a few species are highly dominant [58,60].
Such a large sampling is difficult in Amazonian forests due to logistical challenges. Because the approach used here has been shown to be sufficient to represent the hydraulic trait composition of Amazonian communities [32,43,46], and that the inclusion of more species does not change the overall patterns observed in hydraulic traits [46], we argue that our sampling effort is sufficient to allow us to compare hill and valley communities across our sites.

Hydraulic traits
We used species mean P 50 -the xylem water potential at which 50% of water conductance is lost-as a proxy for tree drought resistance. In Tapajós and PNCV, we measured branch-level P 50 from the relationship between xylem water potential and percentage loss of xylem conductivity (PLC), estimated from percent air discharge (PAD) following the pneumatic method [61]. Data on P 50 for Manaus were obtained from [36,46] which used the same method.
We collected two branches with healthy leaves, longer than 1 m whenever possible, from one to six individuals of each species in the early morning. Branches were bagged in dark, damp plastic bags in the field and then transported to the laboratory, where we performed a second cut at the branch end underwater and left them to rehydrate over a day [32]. We induced embolism by letting the branches dry following the bench dehydration method [62], with progressively longer drying intervals. Using a pressure chamber (PMS 1000; PMS Instruments Co., Albany, OR, USA), we determined the xylem water potential by measuring the leaf water potential, after equilibrating the branch within a black plastic bag for 30 min to 1 h prior to the measurement (following 4). For each water potential measurement, we also measured the air discharged using a Pneumatron device (Pneumatron v1, Plantem-Plant and Environment Technologies, Campinas, Brazil) [63]. We estimated the branch percent loss conductance (PLC) from the PAD obtained from the pneumatic method (ref). To obtain P 50 , we fit a sigmoidal curve to the data, of the form where PAD is percent air discharge, P50 is the water potential Px when 50% of xylem conductance is lost, and Sp is the slope of the curve. Due to the non-linearity of the sigmoidal curve, we first fitted the curve to each branch collected, and discarded branches with spurious PAD values, before aggregating to obtain the individual average, which we later aggregated to species mean.
To gain more insight into the hydraulic strategies in each site, we measured species mean predawn leaf water potential (Ψ pd ), a proxy for rooting depth [64]. Leaf water potentials were measured at the peak of the dry season in July 2019 for the PNCV site and in October 2021 for the Tapajós site. Predawn leaf water potentials were measured with a pressure chamber (PMS 1000; PMS Instruments Co., Albany, OR, USA) from selected individuals, including those sampled for vulnerability to embolism, using two or three healthy, fully expanded leaves. We made Ψ pd measurements between 3 and 6 a.m., always before sunrise.
To place our P 50 observations within a broader context of tree drought resistance, we compiled data from the literature for the tropics (figure 2(a), black dots). The broadest categories contain observations for all Tropical Rainforests (TRR) and Tropical Savannas (TRS), extracted from the Xylem Functional Traits Database [65]. Next, we focused on the two biomes we studied, and compiled data following in topographic classes: Amazônia Hills, Amazônia Valleys and Cerrado Hills (table S4). No data was found on P 50 for Cerrado Valleys. We then compare our measurements ( figure 2(a), red dots) with the compiled dataset. Because most of the data gathered from literature was at the individual level, we also use our data at the individual level (except for Manaus, where only species-level data were available).

Statistical analyses
We first test whether species mean P 50 and Ψ pd are different between hills and valleys at each site (H1). We use the non-parametric Mann-Whitney U test to test whether valleys have greater (less negative) P 50 than hills with a significance level of 0.05. To characterize differential water access, we also tested if valleys have more access to water (valley Ψ pd greater than hill) in Tapajós and PNCV, since data from Manaus was not available. We also tested how our new P 50 measurements compare to the previous values reported in the literature. Tests were conducted using the 'mannwhitneyu' function from the 'scipy' package [66] in Python.
To test H1 and H2, we restrict our analyses to P 50 , and test if average P 50 and its variability differed between valleys and hills and among sites with different rainfall regimes. We fitted five linear models with generalized least squares using site (a proxy for rainfall and seasonality) and topographic position (a proxy for groundwater access) as predictors. Because all water ultimately comes from rainfall, and climatic variability can influence groundwater levels on longer timescales [28], we also include an interaction term between topography and site to test for this combined effect on vegetation.
The simplest model (model 1) was an ordinary linear model that can be written as, where α is the intercept, ε i are the residuals and σ 2 is the standard deviation, which is assumed to be normally distributed and homogeneous among topography and site classes. Model 1 tests for differences in mean P 50 values of species occurring in different valleys and hills in the three sites assuming that each species is independent from each other.
To test whether the variability of P 50 differs between topography and sites, we fitted another three models (2 through 4), which have the same explanatory variables as model 1 but a different error term, where a standard deviation σ j is estimated for the j topographic class (model 2), the j site (model 3) or the j combinations of topography and site (model 4). Thus, each topography class, site or combination of topography and site is allowed to have a different variance, which represents the variability of P 50 in each class. Lastly, phylogenetically close species tend to have similar functional traits because of niche conservatism [67][68][69]. To account for phylogenetic associations of species, we fitted another three models (5 through 7), in which the error structure can be written as, where V is the covariance matrix containing estimates of shared evolutionary history between pairs of species. We build the covariance matrix assuming Brownian motion using the phylogenetic distance between species pairs from a phylogeny with time-calibrated branch lengths generated using the 'V.PhyloMaker' package from R [70]. Thus, models 5, 6 and 7 estimate the variabilities and differences between topographic classes and sites regardless of phylogenetic relatedness, while keeping the error structure associated with site (model 5), topographic position (model 6) or combination of topography and site (model 7). For instance, if a significant difference between hills and valleys in model 1 is no longer significant in model 7, it means that the differences in the P 50 values can be explained by the occurrence of communities from different evolutionary clades in the two topographic classes. On the other hand, if topography is not significant in model 1 but becomes significant in model 7, it means that the P 50 values of the species that occur in the hills and valleys are different from what would be expected simply from the evolutionary distance between the species, indicating possible convergent evolution (homoplasies) between the two environments.
We test the improvement of model fit with increasing complexity using a likelihood ratio test (LRT) and comparing the Akaike information criterion (AIC) values of models with similar numbers of parameters. We fit all models with the 'gls' function of the 'nlme' package in R [71].
To quantify the dissimilarity between the P 50 distributions in each topography-climate class (H2), we calculate the geometric overlap between the P 50 probability density functions of two different classes [72]. Maximum dissimilarity (equal to 1) occurs when there is no overlap between the two curves, and minimum dissimilarity (equal to 0) occurs when there is perfect overlap between both curves. Geometricbased dissimilarities are argued to be more robust than other available methods such as average-based ones, because (i) they do not require standardization of values for comparison with other traits and (ii) are independent of the considered species pool [72,73]. We used Gaussian kernels to fit probability density curves to the species mean P 50 values in each class using the 'gaussian_kde' function of the 'scipy' package in Python [66]. P 50 was significantly less negative in valleys than on hills at the three sites, as stated in our H1 (table S5,  figure 3), and mean P 50 decreases with increasing seasonality for both hills and valleys, indicating higher embolism resistance (table S5, figure 3). Valley Ψ pd was significantly less negative than on hills in Tapajós and PNCV, indicating more access to water for valley species ( figure S1, table 1).

Results
The species measured in our study presented more negative P 50 (greater resistance to embolism) than what was previously reported in the literature for the same biomes and/or topograph  Table 1. Mean and standard deviation (SD) for the water potential associated with 50% conductance loss (P50) and predawn leaf water potential (Ψ pd ) for hills and valleys in our three sites: Manaus, Tapajós and Chapada dos Veadeiros National Park (PNCV). Results of the Mann-Kendall U test are also displayed: asterisks define significant differences (p-value < 0.05) for the test (lesser/greater) and the associated p-value in the Test and p columns.

P50 (MPa)
Ψ pd (MPa) To quantify the effect of rainfall and hydrologic position on P 50 values and their variability (H1 and H2), our fit of Model 1 (simplest ordinary linear model) showed significant differences both between hills and valleys and between sites (table 2, aligned with our H1 and H2). Models 2 and 4 improved the fit when compared to model 1 (LRT: χ 2 = 13.5, d.f. = 1, P < 0.001 and χ 2 = 19.3, d.f. = 5, P = 0.0017 for Model 2 and 4, respectively), with Model 2 having the lowest AIC value of the models without phylogenetic relations (table 3). Model 2 estimated different variances for hills and valleys, which was two times greater in the former. When species' phylogenetic relatedness was considered, the best model estimated different variances for each combination of topographic classes and sites (model 7 in table 2). In this model, both topography and site remained significant predictors. The model suggested that differences are not explained only by species relatedness within communities; however, the mean difference in P 50 between topographic classes decreased while the difference among sites increased in Model 7 when compared to Model 1 (table 2, LRT: χ 2 = 204.2, d.f. = 5, P < 0.0001). Model 7 further corroborated our H2 by showing that valleys have lower P 50 variability than hills (table S5) independent of site. Furthermore, the interaction between topography and sites became significant in Model 7 (table 2). This was due to the reduced predicted difference between hill and valleys as water limitations decreased when accounting for the phylogenetic relatedness (table S5).
Hill and valley P 50 dissimilarity increases as seasonality increases (figure 5), from 0.36 in Manaus to 0.78 in PNCV, aligned with H2 that rainfall totals and seasonality modulate the difference between the P 50 composition in valleys and hills. The dissimilarity also increases between valleys along the rainfall gradient (0.33 between Manaus and Tapajós; 0.36 between Tapajós and PNCV; and 0.56 between Manaus and PNCV valleys). For hills, while the dissimilarity between Manaus and the other sites increases from 0.36 in Tapajós to 0.44 in PNCV, we see high similarity between Tapajós and PNCV hills (0.15), which also have a non-significant difference according to the Mann-Whitney U test performed earlier ( figure 5).

Discussion
Our study demonstrates that rainfall and hydrology interact to modulate species' resistance to embolism at the landscape scale (tables 1 and 2). Reduced rainfall and increased seasonality (figure 2) increase xylem tension through increased leaf transpiration and soil evaporation [6], while hydrology controls differential groundwater access, which can offset the negative effects of atmospheric drought. While topographic influence on P 50 has been documented for specific sites before [36][37][38]50], we now show that this is a significant ecological driver of tree embolism resistance across different rainfall regimes in South America.
Rainfall seasonality increases the dissimilarity in embolism resistance between valleys and hills in accordance with our hypothesis H2 (figure 5). This is driven both by a change in P 50 values (table 1) and variability in each position (table S5), with hills becoming more embolism-resistant (more negative P 50 ) and more variable than their valley counterparts as seasonality increases ( figure 3 and tables 1, S5). While more embolism-resistant species are added to the community (figure 3) with increasing seasonality, they still retain rather vulnerable species. Droughtavoidance strategies such as deep rooting [29], which allows species to access water in the deep soil and reduce investments in xylem resistance [32,74,75], can facilitate the survival of these vulnerable species even in dry climates. The presence of both shallow (more negative Ψ pd ) and deep (less negative Ψ pd ) rooted species in the PNCV (figure S1) highlights the coexistence of different rooting strategies in this drought-stressed area, which we hypothesize can be a mechanism allowing the survival of these vulnerable species in such a stressful environment.
Valley species are consistently less embolismresistant than hills (figure 3, table 1) and have lower One-tailed Mann-Kendall U test results between groups, testing whether P50 values for one group are greater or lesser than another. Note that the comparisons are not bidirectional, and represent rows compared to columns (e.g. Tropical Savannas more resistant than Tropical Rainforests, first row and first column). Blue colors represent a significant, less embolism-resistant result (less negative P50), while red represents a significant, more embolism-resistant result (more negative P50) for the row-column group comparison. Grey colors are non-significant results.   (figure 3,  table S5), while average P 50 values get more negative (higher embolism resistance) when compared to valleys in less seasonal sites. The increased embolism vulnerability in valleys might be associated with the reported tradeoff in which plants maximize their efficiency in water transport from roots to leaves in places of high water availability, to the detriment of xylem resistance to drought [76][77][78][79]. With reduced MAP and longer dry periods, increased atmospheric drought stress might force the most vulnerable trees in valleys beyond progressively more negative xylem pressure thresholds [6]. This would limit the survival of highly vulnerable trees (as seen in Manaus) in more seasonal locations by exposing them to hydraulic failure and carbon starvation [5,[80][81][82]. Moreover, seasonal waterlogging in valleys also imposes a strong abiotic stress and limits the establishment of trees without adaptations to root zone anoxia [83][84][85], which are usually more drought-resistant [84,86].
Our results show that evolutionary links also contribute to the patterns we observed. When phylogenetic relationships are considered (Model 7), both hills and valleys in Manaus are shown to have a larger P 50 range than their counterparts in other sites. The effect of rainfall on P 50 also increased (table 2), indicating convergent evolution acting at the regional (climatic) scale, increasing the difference among the sites from what would be predicted by evolutionary relationships alone. On the other hand, while remaining significant, the mean difference between hills and valleys decreased (decreased F-value of 'topography' , table 2), indicating that communities with species from different lineages are selected on the topographic scale. Moreover, in wetter sites, the phylogenetic structure of hill and valley communities is more distinct than in drier sites. This is evidenced from the stronger predicted difference between hills and valleys at drier sites when considering the phylogenetic relationships among species. These results indicate that environmental filters operate at multiple spatial and temporal scales. On the ecological time scale, the distinct hydrological environments select tree communities with different drought resistance compositions. On the evolutionary time scale, there seems to be a convergent evolution toward lower drought resistance with increasing rainfall [87], at least within the main Angiosperm lineages sampled here.
When compared to reported values in the tropics, our data shows that the variability in P 50 can be much higher than previously thought (e.g. Manaus and Tapajós valleys vs Amazônia valleys, figure 4). Moreover, the median P 50 value of local communities can be highly different from regional assemblages, especially when topography is considered. While part of this difference could be caused by different methods of estimating P 50 , we argue that careful characterization of hydrologic environments and the explicit inclusion of paired hill and valley plots in ecosystem monitoring networks [88,89] are needed to improve the representation of vegetation hydraulic traits heterogeneity across ecosystems.
Overall, our findings highlight the role of topography and rainfall in jointly modulating how evolutionary processes and environmental filtering together shape the resistance of trees to drought [90]. The increased importance of the rainfall-topography interaction in Model 7 (table 2), higher than topographic position alone, suggests that the hydrologic gradient is a stronger environmental filter in more seasonal climates ( figure 5, table 2). In such rainfall regimes, drought stress is high and vegetation productivity, diversity and mortality are strongly driven by water-limitation [43,80,[91][92][93].
We note, however, that our study does not capture the full range of rainfall regimes of South America, but evidence from other sites seem to align with our hypotheses and results. In a high MAP, aseasonal forest in Northwestern Amazônia, no differences were found between hills and valleys for several branchand stomata-level traits [94], although P 50 directly was not measured. Where seasonal water stress is minimal, the relative importance of a shallow WT for drought resistance is likely small, blurring the topographic signal. More research expanding the range of rainfall regimes sampled is urgently needed to improve our understanding of the function of key ecosystems such as Amazônia.

Data availability statements
The data that support the findings of this study are openly available at the following URL/DOI: https:// github.com/caiomattos/erl_embolism.

Funding statement
This research and fieldwork activities were funded by an AGU Horton Research Award to C R C M, a CUAHSI Pathfinder Fellowship to C R C M, a National Geographic Society Early Career Grant