Drought tolerance inferred from functional traits has been shown to explain tree species distributions with respect to rainfall and topography in dry tropical forests, but the extent to which drought tolerance may influence species distributions in wet tropical forests remains unclear.

We measured the leaf water potential at turgor loss point (Ψ_TLP), a leaf hydraulic trait associated with drought tolerance, and representing the leaf water potential threshold for wilting, for 58 abundant tree species in an aseasonal lowland tropical forest in Peninsular Malaysia. We tested whether tree spatial distribution in relation to topographic soil moisture levels (i.e., Topographic Wetness Index; TWI) within a 50-ha forest plot (the ForestGEO Pasoh plot) was related to Ψ_TLP.

We found that 76% of species showed spatial associations with TWI within the 50-ha plot. Species’ spatial association with soil moisture levels was statistically independent of Ψ_TLP.

Even in wet tropical forests defined by no distinct dry season, soil moisture associated with topography exerts a strong influence on tree species distributions, but topographic distributions are not associated with one of the components of species' drought tolerance as quantified with Ψ_TLP.  These findings highlight the sensitivity of tropical forests to variation in soil water levels, regardless of species' capacity to maintain physiological function under drier conditions, leading to a potential shift in community assemblages and structures following extreme events.

Water potential at turgor loss point:

We measured the leaf water potential at the turgor loss point (Ψ_TLP) of the 58 focal tree species. The sampling involved collecting one fully sun-exposed branch from three sun-exposed adult trees (i.e., DBH>25 cm and height above 30m) growing around the field station, located ca. 2 km from the 50-ha plot, between November and December 2018. Indeed, emergent trees are the most vulnerable to drought due to greater hydraulic constraints, evaporative demand, and irradiance (Bennett et al. 2015), while being the most important trees in the forest's water and carbon cycle (Kunert et al. 2017). Cut branches were placed in dark plastic bags with moist paper and transferred immediately to the laboratory at the field station, and prepared for determination of Ψ_TLP via the osmometric method following Bartlett, et al. (2012). The osmometric method allows rapid measurement of the Ψ_TLP, and has been shown to yield values strongly correlated with those obtained from traditional pressure–volume curves in both Amazonian and Malaysian tropical forests (Luiz Ferraresso Conti Junior et al. 2020, Kawai et al. 2025). At the lab, collected branches were re-cut under water to remove possible cavitation at the distal end and placed in buckets with water and covered with dark plastic bags to suppress leaf physiological activity. The branches remained in the buckets for rehydration overnight. On the day following sample collection, we chose two fully developed leaves per branch and cut two 4-mm discs out of each leaf, avoiding major veins. We wrapped the discs in aluminum foil and submerged them in liquid nitrogen. After shock-freezing, the discs were perforated with a dissecting needle 10 to 15 times and immediately placed in the standard 10 μL chamber well of an osmometer (VAPRO 5520, Wescor, Logan, UT, US). We used the auto-repeat mode until an equilibrium was reached (≤ 5 mmol.kg^-1 of variation between measurements). Ψ_TLP was calculated from the osmolality (c_osm) using the calibration equation following Bartlett, Scoffoni, Ardy, et al. (2012):

Then, the mean Ψ_TLP was calculated from the resulting 12 measurement points per species (i.e., 4 measurements per individual and 3 individuals per species). From the species level Ψ_TLP, we calculated the community-weight mean at the quadrat level (cwmΨ_TLP). For each quadrat in the 50ha plot, we summed the products of the species mean measured Ψ_TLP and the species’ proportion of the basal area in the quadrat. Resulting Ψ_TLP indicates the wilting point corresponding to the leaf water potential threshold where turgor and some physiological functions are lost (Bartlett, Scoffoni & Sack 2012): a more negative value maintains leaf turgor under drier conditions.

Characterization of Environmental gradients

1. Topographic Wetness Index (TWI)

We calculated the Topographic Wetness Index (TWI) for the Pasoh 50-ha plot as a proxy of potential soil moisture following the Beven & Kirkby (1979) formulation:

TWI = ln (a / tan(β)), (eq1)

where a is the specific catchment area (i.e., studied cell area, m².m^-1) and β is the local slope angle (radians). TWI quantifies topographic control on hydrological processes and predicts spatial variation in soil moisture availability across the landscape. Indeed, previous studies in dry and wet ecosystems, including Pasoh, showed a strong negative correlation between soil moisture and elevation (Daws et al. 2002, Qiu et al. 2017, Marryanna et al. 2019). The starting digital elevation model (DEM) for the 50-ha plot had a 5 m resolution (Figure S2). However, we estimated TWI using multiple DEM smoothing parameters because topography-derived variables are highly sensitive to the spatial grain of the DEM (Detto et al. 2013, Cushman et al. 2022) and because smoothing helped optimize the balance between noise reduction and the preservation of hydrologically relevant terrain features. We applied a focal smoothing using Gaussian kernels with varying radii (10, 50, and 100 m) and sigma values (5, 50, and 100), testing nine (3 radii x 3 sigma) parameter combinations in total. Gaussian kernels were preferred over circular or rectangular kernels because they preserve terrain gradients through gradual weight decay, avoid artificial flattening at window boundaries, and maintain isotropy in weight distribution—all critical properties for accurate TWI estimation.

For each parameter combination, we: (1) applied focal smoothing on the 5m resolution DEM using normalized Gaussian weight matrices generated with the terra package (version 1.8-54, Hijmans 2025); (2) computed slope in radians from the smoothed DEM using finite-difference methods; (3) filled hydrologic depressions using the Wang & Liu (2006) algorithm implemented in the whitebox package (version 2.4., Lindsay 2016, Wu & Brown 2025) to ensure continuous flow routing; (4) calculated flow accumulation using the D8 deterministic eight-direction flow algorithm, which routes each cell to its steepest downslope neighbor; (5) computed TWI by combining specific catchment area (flow accumulation x cell area) with slope, applying a minimum slope threshold of 1x10^-6 radians to avoid division by zero in flat areas; and (6) average to a grid of 10 x 10 m (Figure 1a). The nine resulting TWI matrices for the 50ha plot were used in the main downstream analyses and for results’ sensitivity assessments. All analyses were carried out in R (version 4.4.2, R Core Team 2024).

 2. Soil habitat classification

We defined soil habitat types based on differences in soil nutrient availability to account for variation in soil chemical properties, widely recognized as a major driver of species distributions across the landscape (John et al. 2007, Zuleta et al. 2020). We selected eight soil nutrients (aluminum [Al], calcium [Ca], iron [Fe], potassium [K], magnesium [Mg], manganese [Mn], sodium [Na], and phosphorus [P]; mg.kg^-1) and pH for habitat classification, as they exhibit high spatial variability and are related to the two primary axes of soil chemical variation in lowland tropical forests (Medina-Vega et al. 2024). Soil samples were collected from the top 10 cm at every other ~ 40 × 40 m quadrat center and, additionally, at distances of 2, 8 or 20 m in a randomly selected direction to capture fine-scale heterogeneity (see details in Medina-Vega et al. 2024). Samples were dried at ambient temperature and sieved at 2 mm. Soil pH was measured using a glass electrode in a 1:2 soil-to-solution ratio in water. Phosphorus (P) was extracted using Bray-1 solution (Bray & Kurtz 1945) and detected by automated molybdate colorimetry on a Lachat QuickChem 8500 (Hach, Loveland, Colorado, USA). The remaining cations were extracted using 0.1 M BaCl₂ (2 hours, 1:30 soil-to-solution ratio) and analyzed by inductively coupled plasma optical emission spectrometry on an Optima 7300 DV (Perkin-Elmer, Shelton, Connecticut, USA; Hendershot et al. 2008).

We interpolated soil chemical data at 10 m resolution using ordinary kriging implemented in the geoR package (v.1.9-4, Ribeiro et al. 2020). Prior to interpolation, raw soil variables were Box-Cox transformed (λ = 0, 0.5, or 1) (John et al. 2007). Assuming isotropy, we fitted a polynomial trend-surface regression of the form: s ~ x + y + x² + y² + x*y, where s is the transformed soil variable and 𝑥, 𝑦 are spatial coordinates (meters). Empirical variograms were calculated from model residuals using the variog function, and covariance parameters (range, sill, and nugget) were estimated by fitting five theoretical variogram models (Gaussian, circular, exponential, spherical, and Cauchy) with the variofit function. The best least-squares model was selected. Kriged predictions were obtained using the krige.conv function, added to the trend-surface predictions, and subsequently back-transformed to produce continuous soil surfaces for habitat classification.

Soil habitat types were identified using spatially constrained hierarchical clustering implemented in the ClustGEO package (version 2.1, Chavent et al. 2021). This approach extends traditional hierarchical clustering by incorporating spatial structure through a mixing parameter (α), which balances environmental dissimilarity against spatial contiguity. The interpolated soil surfaces were converted to 10 m resolution polygons, from which we calculated (1) an environmental dissimilarity matrix based on standardized soil chemical variables, and (2) a spatial distance matrix based on standardized x-y coordinates.

We evaluated α values from 0 (pure environmental constraint) to 1 (pure spatial constraint) in increments of 0.025 for four-cluster solutions. For each α, we computed normalized environmental inertia (quantifying within-cluster chemical heterogeneity) and normalized spatial inertia (quantifying spatial fragmentation). The optimal α was identified as the value that equalized both measurements, representing the optimal trade-off between environmental homogeneity and spatial contiguity; this was determined using a root-finding procedure on the function describing their difference.

Using the optimal α, we performed spatially constrained hierarchical clustering with Ward’s linkage and cut the resulting dendrograms to obtain four habitat classes. The final soil habitat map delineates spatially coherent zones characterized by distinct soil chemical signatures: (1) acidic soils enriched in Mg and Mn, (2) high-pH soils with elevated nutrient concentrations, (3) nutrient-poor soils enriched in Fe, and (4) Ca- and Na-rich soils (Figure 1b, Figure S3).