Machine-learning models of δ¹³C and δ¹⁵N isoscapes in Amazonian wood

Illegal logging is one of the most prevalent environmental infractions in the Amazon, led by organized networks that cause substantial ecological and economic impacts. Official control mechanisms, such as Brazil’s Forest Origin Document (DOF), remain vulnerable to the fraudulent manipulation of virtual timber credits and inconsistencies in digital traceability. These deficiencies highlight the need for independent, scientifically based methodologies for timber traceability that can support law enforcement and ensure reliable provenance verification. Here, we tested whether the isotopic composition of carbon (δ¹³C) and nitrogen (δ¹⁵N) in wood can trace Amazonian timber origin. We developed basin-wide δ¹³C and δ¹⁵N isoscapes using machine-learning models to predict spatial variability. A total of 571 trees from 47 sites were analyzed for both isotopes. Tree disks or wedges were sampled from the basal trunk, sectioned transversely, and sub-sampled from heartwood to near the sapwood boundary to obtain whole-tree isotopic composition. The δ¹³C and, more strongly, the δ¹⁵N values exhibited substantial within-site heterogeneity, indicating individual-level physiological controls, interspecific differences, and/or small-scale environmental variation influencing isotope fractionation. Despite these sources of noise, isotopic values showed independent and predictable spatial patterns across the basin (R²=0.67 for δ¹⁵N and R²=0.60 for δ¹³C). Nitrogen isotopes were primarily controlled by edaphic factors, while carbon isotopes revealed a broad longitudinal gradient linked to climate. Together, these isotopic markers provide complementary information for basin-scale timber provenancing and form a robust, high-resolution framework for Amazon-wide traceability.

1 Introduction

The Amazon rainforest plays a critical role in regulating climate processes on continental and global levels, influencing rainfall patterns from the Atlantic Ocean to the Andes and to south-central Brazil (Salati et al., 1979). However, the region has lost more than one-fifth of its forest cover since the 1980s (INPE, 2023), driven largely by land grabbing, agricultural expansion, and illegal logging (Nepstad et al., 2014; Matricardi et al., 2020; Ferrante et al., 2021; Silva Junior et al., 2020; Saraiva, 2021). Globally, 15 %–30 % of traded timber is estimated to be illegal, reaching 50 %–90 % in tropical regions (INTERPOL, 2019). In the Amazon, roughly 40 % of timber extraction between August 2020 and July 2021 occurred without authorization, mainly on private properties and partly in Indigenous lands and protected areas (Valdiones et al., 2022). Although governance initiatives such as the Non-Prejudicial Extraordinary Report (NDF) under CITES Appendix II aim to improve legality verification (IBAMA, 2024), deficiencies persist, enabling misidentification, volume manipulation, and “legalization” of timber of illicit origin (Brancalion et al., 2018; CNI, 2018).

In Brazil, provenance control based on the Forest Origin Document (DOF) and Forest Guides remains primarily documentary and legally fragile, enabling inconsistencies along the logistics chain and undermining trust in official traceability (Andrade et al., 2023; Franca et al., 2023). Therefore, independent verification methods grounded in material evidence are needed (Brancalion et al., 2018). Several analytical approaches have been developed to verify the geographic origin of timber across both temperate and tropical regions, including dendrochronological and anatomical analyses, multi-elemental chemistry, and genetic fingerprinting (Akhmetzyanov et al., 2019; Delmás et al., 2020; Hornink et al., 2025; Boeschoten et al., 2025). These methods provide complementary lines of evidence and can be combined to strengthen wood traceability frameworks.

Among them, stable isotopes serve as intrinsic markers of provenance because they can imprint distinct isotopic compositions in organic tissues growing at different locations. Spatial variations in temperature, water availability, altitude, soil properties, and lithology drive isotope variations in the environment, and those variations are propagated to plants and ecosystems (Leavitt and Long, 1984; Deleens et al., 1994; McCarroll and Loader, 2004; Fry, 2006; Leavitt and Roden, 2022).

In C₃ plants, which include all tree species, carbon isotopes are assimilated from atmospheric CO₂ that diffuses from the atmosphere to the leaf carboxylation sites, where the rubisco enzyme converts them into carbohydrates. The $p_{\mathrm{i}}/p_{\mathrm{a}}$ ratio is influenced by both CO₂ diffusion into the leaf and carboxylation rate and is directly linked to carbon isotope fractionation (Farquhar et al., 1982, 1989). Despite their identical photosynthetic pathways, C₃ plants can exhibit a broad range of δ¹³C values, particularly in forested ecosystems, influenced by the “canopy effect” and by stomatal conductance (Farquhar et al., 1989; Medina and Minchin, 1980; Ometto et al., 2006). Plants growing in closed-canopy forests tend to show lower (more negative) δ¹³C values because of the incorporation of respired CO₂, which has lower ¹³C content. In contrast, plants in drier environments often display higher (more positive) δ¹³C values, as reduced stomatal conductance limits CO₂ diffusion (Farquhar et al., 1989; Martinelli et al., 2021). However, carbon isotopes have rarely been used as a marker of provenance because intra-site δ¹³C variations are large (Martinelli et al., 1998; Paredes-Villanueva et al., 2022). At a given site, co-occurring tropical forest tree species often show large δ¹³C differences, reflecting variations in successional status, canopy position, and consequently in water use efficiency (WUE) (Guehl et al., 1998; Bonal et al., 2000).

Nitrogen isotopes in plants are primarily controlled by the source of nitrogen assimilated because uptake of nitrogen by plants leads to negligible fractionation (Högberg, 1997). However, soil nitrogen cycling and soil microbial processes strongly fractionate nitrogen isotopes, and those δ¹⁵N variations are propagated to plants (Robinson, 2001; Craine et al., 2015). Among the key processes fractionating nitrogen isotopes, nitrification and the subsequent leaching of soluble nitrates out of soils in wet environments can leave soils with higher δ¹⁵N values (Mariotti et al., 1981; Craine et al., 2015). Denitrification through the diffusion of ¹⁵N-depleted gases also leaves soils with higher δ¹⁵N values (Mariotti et al., 1981; Wang et al., 2018). Ectomycorrhizal fungi also play a key role in controlling δ¹⁵N values in plants by influencing how tightly nitrogen cycles in soils, thereby limiting or amplifying fractionation processes (Craine et al., 2009). Nitrogen-fixing legumes, in principle, can also lower foliar δ¹⁵N values because atmospheric N₂ has a δ¹⁵N close to 0 ‰ (Evans, 2001; Vitousek et al., 2002).

In tropical forests, however, soils and foliage are generally higher in ¹⁵N compared with other biomes, reflecting an open nitrogen cycle dominated by losses (Martinelli et al., 1999). In the Amazon, this pattern is especially evident; although nitrogen-fixing trees are present, in mature terra firme forests they often contribute little to the nitrogen balance because most legumes rely on soil nitrogen rather than atmospheric fixation (Ometto et al., 2006; Nardoto et al., 2008). Moreover, ectomycorrhizal associations are relatively rare, with arbuscular mycorrhizae predominating (Corrales et al., 2018). Since AM fungi exert weaker control on isotopic fractionation, Amazonian plants tend to reflect the δ¹⁵N values of the soil nitrogen pools more directly. δ¹⁵N values in the Amazon are further shaped by edaphic properties and landscape position, which regulate the relative importance of processes such as nitrification, leaching, and denitrification (Martinelli et al., 1999; Amundson et al., 2003; Ometto et al., 2006; Nardoto et al., 2014).

Although both isotopes respond to fine-scale environmental and physiological conditions that may introduce some within-site heterogeneity, coherent large-scale patterns of δ¹³C and δ¹⁵N have been consistently observed across South America, particularly within the Amazon Basin (Martinelli et al., 2021). In the Brazilian Amazon, soil δ¹⁵N values are higher in the east than in the west (Sena-Souza et al., 2020), with lower isotopic compositions occurring in the nutrient-rich, younger soils of western Amazonia and higher values in the highly weathered, nutrient-poor soils of the eastern region (Martinelli et al., 1999; Nardoto et al., 2014). Across the continent, foliar δ¹³C values are lowest in tropical forests, particularly within the Amazon Basin (Powell et al., 2012). Within the basin, δ¹³C values become progressively less negative from the humid core areas toward the drier transition zones and forest peripheries (Ometto et al., 2002; Martinelli et al., 2007). This coherent spatial structuring highlights the strong environmental control on isotopic composition and underscores the potential of combining carbon and nitrogen isotopes as complementary intrinsic tracers for determining provenance within tropical ecosystems.

Isotopic maps, or isoscapes, render these spatial patterns more explicit and are extensively used in the fields of ecology, biogeochemistry, food authentication, and wildlife tracking (Bowen, 2010; Bataille and Bowen, 2012; Watkinson et al., 2020, 2022; Bataille et al., 2021; Reich et al., 2021; Le Corre et al., 2025). In forensic contexts, they can provide indications or material evidence recognized in judicial processes (Ehleringer and Matheson Jr., 2010; Chesson et al., 2018; Meier-Augenstein, 2017). Among the various methods used to construct isoscapes, recent studies have increasingly employed machine-learning techniques to model complex and nonlinear isotopic patterns across landscapes (Barberena et al., 2021; Reich et al., 2024; Holt et al., 2025).

The absence of carbon and nitrogen isoscapes for Amazonian wood precludes the potential use of these isotopes for origin attribution methodologies, along with other isotopes (e.g. strontium, sulfur and oxygen) or other provenance markers (e.g. genetics). We aim at testing whether δ¹³C and δ¹⁵N in wood across the Amazon are predictable and contain provenance information for Amazonian wood at the basin scale, and whether taxonomic effects are overridden by environmental and biogeochemical controls at this scale. We first generated a dataset of carbon and nitrogen isotopes in wood across the Amazon. We then leveraged machine-learning approaches to predict the spatial variability of carbon and nitrogen isotopes across the landscape.

2 Materials and methods

2.1 Study area and sampled species

From a representative collection of 571 trees encompassing 75 genera and 25 botanical families distributed across 47 sites throughout the Amazon, taxonomic identification was performed at the genus level because species-level classification is particularly challenging in tropical forests with high floristic diversity. A subset of samples originated from naturally fallen timber found in the field, which limited the possibility of taxonomic identification. These were classified as “unknown”, representing approximately 5 % of the total dataset (29 individuals).

The selection of species prioritized taxa with high commercial value in the Brazilian Amazon, specifically those identified as most exploited by IMAFLORA (Andrade et al., 2022). Due to logistical limitations in the sampling of certain target species, additional genera with lower economic significance were included based on availability through partner research groups. Nonetheless, the majority of the sampled individuals belong to genera acknowledged for their commercial importance in the region (Fig. S1 in the Supplement).

The selection of the 47 sampling sites aimed at maximizing representation of the region's spatial and ecological variability, with sites encompassing a broad geographic distribution across the Brazilian Amazon (Fig. 1). In practice, while some samples were precisely georeferenced, the majority of the samples were collected within a site to simplify collection for our forestry partners. In general, sites were defined as zones of a few square kilometers where sampling of individual trees was conducted along the road within a few kilometers (min =0 km, max =215 km, median =0.9 km).

Figure 1 Sampling sites across the Brazilian Amazon. (a) South America, with emphasis on the Brazilian Amazon region. (b) Distribution of sampled individuals across 47 sites in the Brazilian Amazon. Numbers correspond to site identifiers listed in Table 1. The map highlights the states where sampling took place, with grey lines indicating state boundaries. Yellow symbols represent sampled sites, while blue symbols (sites 20, Feijó, and 24, Itapuã do Oeste) indicate the locations analyzed in detail in this study. Symbol size is scaled according to the number of individuals sampled per site (see legend). South America boundaries were obtained from Natural Earth (2023), and the Amazon biome polygon from Assis et al. (2019).

2.2 Sample collection and preparation

The minimum diameter of sampled trees adhered to the legal harvesting threshold (DBH ≥ 50 cm) as stipulated by Brazilian regulations (IBAMA, 2006; CONAMA, 2009). To ensure the geographic traceability of the samples, the majority of sites were situated within protected areas, while the other sites comprised certified sustainable forest management areas and long-term forest research plots. Geographic coordinates of the sampling sites and individual trees are provided in Table S3, which is available in the OSF repository linked to this paper instead of being included directly in the Supplement.

From each tree, a cross-sectional disk or wedge was excised at the basal part of the trunk, preserving the complete radial profile inclusive of the pith, heartwood, and sapwood. The disks measured approximately 5–6 cm in thickness. In the laboratory, each specimen was transversely sectioned into a strip with dimensions of 1.5 cm in width and 2.5 mm in thickness, and subsequently oven-dried at 70 °C for 24 h, following standard protocols for tropical wood anatomical and densitometric preparation (Quintilhan et al., 2021). Subsequent sampling was guided by the work of Batista et al. (2025), who evaluated isotopic variation along the radial axis by excising five 1.5 × 2 cm segments at distinct anatomical positions (pith, heartwood, heartwood–sapwood transition, and sapwood). Their results showed that the heartwood–sapwood transition provides a reliable approximation of mean isotopic composition across the entire radius. Based on this evidence, we targeted this specific segment for the present study (Fig. 2), as it combines analytical representativeness with practical advantages for forensic applications, since it is typically preserved in processed timber, whereas the sapwood is often removed during industrial processing.

Figure 2 Field sampling procedure. (1) Cross-sectional disk excised at the base of the trunk, preserving the anatomical structure of the wood. (2) Disk sectioned with approximately 5–6 cm thickness. (3) Transversal sectioning of the wood strip from the disk and subsequent sampling of the transition zone between heartwood and sapwood along the radial axis for isotope analysis.

2.3 Isotope analyses

For the purpose of conducting isotopic analysis, the samples were meticulously ground using a RETSCH ZM 300 mixer mill to ensure thorough homogenization. All procedures were executed within a climate-controlled environment to minimize the potential for contamination. The resultant powdered samples were subsequently weighed using an analytical balance, with masses spanning from 2.7 to 3.0 mg to ascertain the requisite quantity for analysis. Thereafter, the pulverized material was encapsulated within tin capsules, which are appropriate for the determination of ¹³C $/$ ¹²C and ¹⁵N $/$ ¹⁴N isotopic ratios.

The isotopic analyses were performed utilizing a Carlo Erba 1110 elemental analyzer (Milan, Italy), coupled a isotope ratio mass spectrometer (Delta plus 2, Thermo Fisher Scientific, Bremen, Germany) for the determination of isotopic ratios. Analyses were undertaken at the Isotopic Ecology Laboratory of the Center for Nuclear Energy in Agriculture (LEI/CENA/USP, Brazil).

The isotopic ratio (R) was defined as the ratio of the rare to the abundant isotope in the sample, expressed relative to an international reference standard. Isotopic values were calibrated against the Vienna Pee Dee Belemnite (VPDB) and AIR scales for carbon and nitrogen, respectively. For carbon and nitrogen isotope analyses, internationally certified reference materials NBS-22 (${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=-\mathrm{30.03}\pm \mathrm{0.04}\mathrm{‰}$ VPDB), IAEA-N-1 (${\mathit{\delta }}^{\mathrm{15}}\mathrm{N}=+\mathrm{0.43}\pm \mathrm{0.07}\mathrm{‰}$ AIR), and IAEA-N-2 (${\mathit{\delta }}^{\mathrm{15}}\mathrm{N}=+\mathrm{20.41}\pm \mathrm{0.12}\mathrm{‰}$ AIR) were employed for calibration. The internal laboratory standard consisted of a sugarcane material routinely analyzed as a blind quality-control sample. This internal standard was measured repeatedly throughout each analytical sequence to monitor instrument stability and analytical precision. Reported values remained within the expected range, with standard uncertainties of ±0.10 ‰ for δ¹³C and ±0.15 ‰ for δ¹⁵N.

Isotopic composition is expressed in delta (δ) notation, with values reported in per mil (‰) relative to international reference standards (e.g., VPDB for carbon, AIR for nitrogen), following IUPAC recommendations (Coplen, 2011; Skrzypek et al., 2022), and defined as:

$$\begin{array}{}\text{(1)}& {\mathit{\delta }}^{i/j}E={\displaystyle \frac{R_{\text{sample}}}{R_{\text{standard}}}}-\mathrm{1}\end{array}$$

where R represents the ratio of the heavy (iE) to the light (jE) isotope of element E (e.g., ${}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ or ${}^{\mathrm{15}}\mathrm{N}{/}^{\mathrm{14}}\mathrm{N}$).

2.4 Intra-site variability

All data processing, statistical analyses, and figure generation were performed in R version 4.3.3 (R Core Team, 2025). While the goal of this study is to primarily assess the inter-site variability and the potential of using carbon and nitrogen isotopes for tracing wood origin, we first examined intra-site isotopic variance to evaluate the portion of the total variance that is not dependent on spatial autocorrelation. In this study, we sampled trees of different ages, species, heights, and rooting depths, which are known to influence carbon and nitrogen isotopes as shown by Batista et al. (2025). Additionally, the loose definition of “site” in this study, as a zone of several square kilometers potentially encompassing distinct geomorphology, geology, microclimates, and environments, would inherently lead to some intra-site variability.

For each site, we calculated the median, mean, standard deviation (SD), and root mean square error (RMSE) for δ¹³C and δ¹⁵N values. Boxplots were generated to illustrate the distribution of values within each site. These descriptive statistics provide the basis to evaluate the intra-site uncertainty for each isotope and the spatial resolution limitations of their respective predictive models.

To further disentangle geographic and taxonomic sources of variability, we conducted specific analyses on three genera (Cedrela, Handroanthus, and Dipteryx) that co-occurred in two well-sampled sites (Feijó and Itapuã do Oeste). Site effects were tested using independent-sample t-tests, while genus effects within sites were assessed using one-way ANOVA followed by pairwise comparisons with false discovery rate (FDR) correction. We also calculated 95 % confidence intervals (CI) around mean values to quantify within-site dispersion at the genus level. These combined analyses allowed us to evaluate whether isotopic variability was more strongly influenced by taxonomic identity or geographic origin.

2.5 Inter-site variability and isoscape development

In order to predict δ¹³C and δ¹⁵N values in Amazonian wood, a multivariate regression approach was employed, following methodologies established in prior isoscape studies (e.g., Bataille et al., 2018, 2020; Sena-Souza et al., 2020; Reich et al., 2021; Le Corre et al., 2025; Martinelli et al., 2025). The process involved four main steps: (1) collection of auxiliary variables; (2) application of a Random Forest regression algorithm; (3) assessment of model performance through cross-validation and prediction error metrics; and (4) spatial application of the model to produce predictive isoscapes and corresponding uncertainty maps. The subsequent sections elaborate on each step in detail.

2.5.1 Auxiliary variables

To predict the spatial variability of δ¹³C and δ¹⁵N values within Amazonian wood, we assembled an extensive array of 74 spatially explicit covariates that encompass critical environmental, climatic, edaphic, geophysical, and biological factors (Table S1 in the Supplement). These variables were selected based on their established direct or indirect associations with plant isotopic composition and their ability to represent processes influencing carbon and nitrogen assimilation within tropical forest ecosystems.

The dataset includes soil physicochemical properties such as pH, cation exchange capacity, clay content, bulk density, total nitrogen, and organic carbon stocks, which govern nitrogen mineralization, microbial activity, and root nutrient uptake factors intrinsically linked to δ¹⁵N variation (Amundson et al., 2003; Nardoto et al., 2008; Craine et al., 2015; Savard and Daux, 2020; Brunello et al., 2024). Climatic variables, including temperature, precipitation, vapor pressure, aridity index, and potential evapotranspiration, affect stomatal conductance and photosynthetic discrimination, which are primary drivers of δ¹³C variability in C₃ plants (Farquhar et al., 1982; Locosselli et al., 2013; Van der Sleen et al., 2014; Martinelli et al., 2021).

We also incorporated remotely sensed indicators of vegetation structure and function, such as canopy greenness, leaf area, net primary productivity, and carbon assimilation. These proxies reflect photosynthetic capacity and ecosystem productivity, consequently impacting plant nutrient demand and isotopic signatures (Ometto et al., 2006; Cernusak et al., 2013). Wood density, an ecophysiological trait related to water-use efficiency, was included as a proxy for genus-specific physiological strategies (Halder et al., 2024).

Furthermore, we examined topographic attributes such as elevation, geographic gradients like proximity to the coastline, and an array of geological and geophysical variables encompassing bedrock age, porosity, and the atmospheric deposition of dust and sea salt. A selection of raster layers utilized in this investigation was acquired in a pre-processed format from a geospatial database assembled by Bataille et al. (2018, 2020, 2021), originally designed for global and regional isoscape modeling of strontium isotopes. Additional covariates were independently acquired and processed from publicly accessible datasets to ensure comprehensive spatial coverage and thematic complementarity. All environmental geospatial products were resampled and reprojected to the equivalent WGS84–Eckert IV projection, with a spatial resolution of 1 km², to standardize the raster layers and mitigate area distortions at a global scale.

2.5.2 Random forest regression and geospatial predictions

Random Forest constitutes a tree-based machine learning algorithm that constructs multiple decision trees by employing bootstrap sampling alongside random feature selection (Breiman, 2001). It synthesizes their outputs to estimate a response variable while obviating the necessity for data transformation or presumptions regarding distribution or residual variance (Bataille et al., 2020). In accordance with the general Random Forest regression methodology established by Bataille et al. (2018) for ⁸⁷Sr/⁸⁶Sr isoscapes and implemented using the caret package (Kuhn, 2008), we refined and applied this approach to model the spatial distribution of δ¹³C and δ¹⁵N.

Before the initiation of model training, environmental variables were extracted from raster layers corresponding to each individual sampling point and submitted in full to the VSURF (Variable Selection Using Random Forest) algorithm. After the selection process, a Pearson correlation analysis was performed among the selected variables to identify highly correlated pairs (r>0.90), which were considered redundant. Additionally, some variables with lower predictive contribution were excluded to improve model parsimony and reduce error. Although Random Forest algorithms are inherently robust to multicollinearity, these post-selection refinements aimed to enhance both the ecological interpretability and the predictive performance of the final model.

Discrepancies in spatial georeferencing were noted among the sampling sites; in some instances, individual trees were assigned unique coordinates, while in other cases, all individuals from a single site shared a common coordinate. To mitigate spatial pseudoreplication, data were aggregated at the site level. For each site, median values of latitude, longitude, δ¹³C, and δ¹⁵N were computed, resulting in 47 unique site-level entries.

Subsequently, site-level entries were utilized to derive the values of 74 environmental variables from raster datasets, culminating in a comprehensive regression matrix. To discern the most pertinent predictors of δ¹³C and δ¹⁵N, the VSURF algorithm (Genuer et al., 2015) was employed, executing 3000 trees at each phase. Variable selection was predicated on the “prediction” phase of the algorithm, which isolates the subset of predictors that most effectively reduce prediction error.

The designated predictors were subsequently employed to construct Random Forest regression models using the caret package, which incorporated hyperparameter tuning and employed repeated 10-fold cross-validation (with five repetitions), with 80 % of the data allocated for training during each iteration. Model efficacy was assessed using the coefficient of determination (R²), root mean square error (RMSE), and mean absolute error (MAE).

To elucidate the contribution of each variable to the model predictions, we calculated variable importance metrics grounded in the framework of the Random Forest algorithm. Subsequently, for all selected covariates, partial dependence plots (PDPs) were generated to illustrate the marginal effects of each geospatial predictor on the isotopic ratios of δ¹³C and δ¹⁵N (Friedman, 2001). Additionally, to investigate potential interaction effects between environmental variables on isotopic values, bivariate PDPs were constructed using the two most influential covariates. These plots facilitate a more comprehensive interpretation of the combined influence of predictors on the predicted isotopic values, thereby enhancing the ecological understanding of the observed spatial patterns.

Utilizing the refined Random Forest model and the selected predictors, isoscapes were developed to depict the mean predicted δ¹³C and δ¹⁵N values in wood across the Amazon region. To ensure spatial consistency with regional datasets, the final maps were reprojected to the SIRGAS 2000 geographic coordinate system (EPSG:4674), which is suitable for South America. To appraise the spatial uncertainty associated with these predictions, a crucial process in evaluating the reliability of isoscapes, we employed the Quantile Regression Forest (QRF) model, as delineated by Meinshausen (2006) and Le Corre et al. (2025). This algorithm facilitates the estimation of prediction intervals based on the distribution of outputs from the ensemble of decision trees. For each pixel, the 15.9th and 84.1st quantiles were extracted, roughly corresponding to ±1 standard deviation under the assumption of normally distributed residuals. Due to computational memory constraints, predictions were performed in successive blocks of 100 000 pixels containing complete environmental data. Spatial uncertainty was subsequently computed as half the difference between the upper and lower quantiles, resulting in a continuous raster that represents the standard error of the prediction for each pixel. This methodology enabled the development of high-resolution uncertainty maps (1 km²), which can be integrated into probabilistic geographic assignment models, thereby augmenting the robustness and reliability of forensic and ecological applications based on δ¹³C and δ¹⁵N isotopes.

3 Results

3.1 Inter- and intra-site variability

The isotopic composition of all 571 wood samples adhered to a normal distribution for both δ¹³C and δ¹⁵N (Fig. S2). Mean δ¹³C was $-\mathrm{28.1}\pm \mathrm{1.4}$ ‰, with a 95 % confidence interval ranging from −28.2 ‰ to −28.0 ‰, spanning a total range of 9.7 ‰ (−33.7 ‰ to −24.0 ‰). Conversely, δ¹⁵N averaged 4.0±2.3 ‰, with a 95 % confidence interval of 3.8 ‰ to 4.2 ‰, and varied from −2.7 ‰ to 11.2 ‰, resulting in a total amplitude of 13.9 ‰.

Site-level variability in δ¹³C (Fig. 3a) and δ¹⁵N (Fig. 3b) is summarized in Table 1. For δ¹³C, mean site values were lowest in Juruá ($-\mathrm{30.5}\pm \mathrm{1.4}$ ‰) and Atalaia do Norte ($-\mathrm{30.4}\pm \mathrm{1.1}$ ‰), both located in western Amazonia, and highest in Rurópolis ($-\mathrm{26.7}\pm \mathrm{1.0}$ ‰) and Ferreira Gomes ($-\mathrm{26.6}\pm \mathrm{0.9}$ ‰), in the eastern and northeastern basin. The greatest intra-site δ¹³C variability occurred in Alvorada d’Oeste (SD ±1.9 ‰; RMSE = 1.63 ‰; n=4), Uruará (SD ±1.7 ‰; RMSE = 1.53 ‰; n=6), and Pauini (SD ±1.5 ‰; RMSE = 1.50 ‰; n=41), whereas Porto Esperidião (SD ±0.3 ‰; RMSE = 0.28 ‰; n=5), Santa Maria das Barreiras (SD ±0.4 ‰; RMSE = 0.35 ‰; n=9), and Cáceres (SD ±0.4 ‰; RMSE = 0.38 ‰; n=12) exhibited the lowest dispersion.

For δ¹⁵N, mean site values varied widely, with the lowest means recorded in Comodoro (0.6±2.5 ‰) and Feijó (1.6±1.6 ‰), and the highest in Candeias do Jamari (7.5±1.9 ‰). The highest individual value (9.0 ‰) was observed in Monte Negro and was based on a single sampled tree. The greatest intra-site variability occurred in Vitória do Xingu (SD ±3.0 ‰; RMSE = 2.86 ‰; n=9), Porto Velho (SD ±2.9 ‰; RMSE = 2.71 ‰; n=12), and Nova Brasilândia d’Oeste (SD ±2.8 ‰; RMSE = 2.48 ‰; n=6), whereas Pimenta Bueno (SD ±0.0 ‰; RMSE = 0.01 ‰; n=3), Castanheiras (SD ±0.5 ‰; RMSE = 0.36 ‰; n=3), and Alvorada d’Oeste (SD ±0.5 ‰; RMSE = 0.37 ‰; n=4) displayed minimal dispersion.

Table 1 Site-level summary of δ¹³C and δ¹⁵N in wood samples. Includes site number (n), site name, state, geographic coordinates, number of individuals, and the median, mean ± SD, and RMSE for both isotopes.

Note: “–” indicates sites without isotopic data;“NA” marks unavailable statistics due to single observations.

Figure 3 Site-level distributions of δ¹³C and δ¹⁵N in wood, ordered by site median. Numbers along the y-axis correspond to site names listed in Table 1. Gray dots show individual trees; vertical lines indicate medians; boxes show interquartile ranges (IQR); whiskers extend to 1.5× IQR. Sites were grouped by (a) mean annual relative humidity (≥70 % vs. <70 %), based on the humidity threshold identified in this study, and by (b) mean geological age (<539 Ma vs. ≥539 Ma; Cohen et al., 2025).

Sample size across sites ranged from 1 to 50 individuals. In general, sites with larger sample sizes tended to exhibit greater within-site variability, likely reflecting increased ecological or taxonomic heterogeneity. For instance, Pauini (n=41) and ZF2 (n=17) presented relatively high δ¹³C standard deviations (both ±1.5 ‰), while Porto Velho (n=12) exhibited elevated δ¹⁵N variability (±2.9 ‰). However, this pattern was not universal. Comodoro (n=4) showed high δ¹⁵N variability (±2.5 ‰), whereas Feijó (n=30) displayed comparatively low intra-site variation (±1.6 ‰). Notably, the site with the largest sample size, Itapuã do Oeste (n=50), exhibited only moderate within-site variation for both isotopes (δ¹³C SD ±1.0 ‰; RMSE = 0.97 ‰; δ¹⁵N SD ±1.3 ‰; RMSE = 1.29 ‰), indicating that sample size alone does not explain intra-site variability.

To further evaluate intra- and inter-generic variability, we focused on three genera (Cedrela, Handroanthus, and Dipteryx) that co-occurred in two well-sampled sites, Feijó and Itapuã do Oeste. For δ¹³C, differences among genera were more pronounced than differences between sites (Feijó: $-\mathrm{27.76}\pm \mathrm{1.54}$ ‰; Itapuã do Oeste: $-\mathrm{27.76}\pm \mathrm{1.24}$ ‰; Fig. 4a), with no significant site effect (t=0.016, df = 29, p=0.987). In Feijó, all three genera differed significantly (p<0.02), whereas in Itapuã do Oeste only Cedrela and Handroanthus contrasted significantly (p=0.0114). Handroanthus consistently showed the highest δ¹³C values, Cedrela the lowest, and Dipteryx intermediate, with substantial within-site dispersion (Table S2).

In contrast, for δ¹⁵N the site effect was more pronounced than the taxonomic effect (Fig. 4b). Individuals from Itapuã do Oeste exhibited systematically higher δ¹⁵N values compared to those from Feijó (5.67±1.19 ‰ vs. 2.25±1.67 ‰), a highly significant difference ($t=-\mathrm{6.30}$, df = 29, p<0.0001). This site-driven pattern was consistent across all three genera, with Handroanthus, Cedrela, and Dipteryx each showing substantially higher δ¹⁵N values in Itapuã do Oeste than in Feijó (Table S2). Within Feijó, inter-generic contrasts were evident, with Handroanthus differing significantly from both Cedrela and Dipteryx, whereas in Itapuã do Oeste no significant inter-generic differences were detected (all p>0.5). Within-site dispersion at the genus level remained variable (Table S2).

Figure 4 δ¹³C and δ¹⁵N values in wood by genus (Handroanthus, Cedrela, and Dipteryx) across two sites (Feijó and Itapuã do Oeste). (a) δ¹³C and (b) δ¹⁵N distributions. Gray dots represent individual trees; vertical bars indicate medians; boxes represent interquartile ranges (IQR); whiskers extend to 1.5× IQR. Different lowercase letters indicate significant differences among genera within the same site according to Tukey's HSD test (p<0.05, FDR-adjusted).

3.2 Random Forest regression and geospatial predictions

3.2.1 δ¹³C isoscape

After variable selection using the VSURF algorithm, two climate-related covariates associated with atmospheric moisture were identified as the main predictors of δ¹³C in wood: mean annual vapor pressure (r.vapor) and mean annual relative humidity (r.hurs_mean), both showing similar importance in node purity within the model. Following n-fold cross-validation (see Materials and Methods, Sect. 2.5.2), the Random Forest model explained 60 % of the variance in observed δ¹³C values, with an RMSE of 0.59 ‰ and a MAE of 0.48 ‰.

The scatter plot of observed versus predicted δ¹³C values reveals an approximately linear relationship (Fig. 5a). However, the fitted regression line (black) deviates from the identity line (dashed), indicating a tendency to overestimate δ¹³C at the lower end of the distribution (more negative δ¹³C values) and to underestimate at the upper end. This systematic deviation can be referred to as a model bias, characterized by regression toward the mean. The highest density of points occurs within the intermediate δ¹³C range (−28 ‰ to −27 ‰), with fewer observations toward both extremes.

Spatially, residuals are heterogeneously distributed, with both positive and negative values scattered across the basin (Fig. S3). This lack of a geographic pattern suggests that the model does not exhibit systematic spatial bias and performs well at a regional scale, although local environmental factors not captured by the predictors may still influence δ¹³C in wood.

Partial dependence plots (PDPs) indicate nonlinear relationships between δ¹³C and the two selected predictors (Fig. 5b). For relative humidity, δ¹³C values remain stable below ∼63 %, increase between ∼64 %–68 %, and decline above ∼69 %, indicating more negative δ¹³C under high-humidity conditions. A similar pattern occurs for vapor pressure: δ¹³C is lower below ∼2.6 kPa, increases between ∼2.6–2.8 kPa, and decreases again beyond ∼2.85 kPa. The interaction plot (Fig. 5c) shows that the highest predicted δ¹³C values occur under combined conditions of vapor pressure between ∼2.7–2.85 kPa and relative humidity between ∼63 %–68 %, with declines outside this range. The atmospheric vapor pressure and relative humidity layers are provided in the Supplement for interpretation (Fig. S4).

Based on the final model trained with 47 sites, a δ¹³C isoscape was generated for the Amazon region (Fig. 6a). Predicted values range from −32 ‰ to −25 ‰, with lower values concentrated across much of Amazonas and along the Amazon River toward the northeast, and higher values predominantly in the southern, southeastern, and parts of the central-eastern basin. This spatial structure reflects climatic gradients, particularly atmospheric vapor pressure and relative humidity. Predictive uncertainty (Fig. 6b) is higher in areas with sparse sampling, especially in the eastern, central, and northwestern Amazon, and lower in regions with denser sampling coverage, particularly in the south.

Figure 5 Random Foest regression performance and climatic partial dependence plots (PDPs) for δ¹³C in wood. (a) Relationship between observed and predicted δ¹³C values, with a fitted regression line (black) and the 1:1 reference line (dotted). (b) PDPs of mean annual relative humidity (left) and mean annual vapor pressure (right) on predicted δ¹³C. Solid black lines represent the average predicted response when varying each predictor independently. Histograms below each plot indicate the distribution of observed values for the corresponding variable, providing context on sampling density. (c) Interaction surface between relative humidity and vapor pressure showing their combined influence on predicted δ¹³C. The color gradient represents predicted δ¹³C values (‰), with warmer tones indicating higher values and cooler tones indicating lower values. Black contour lines delineate zones of similar predicted values.

Figure 6 δ¹³C isoscape and associated uncertainty. (a) Predictive isoscape of δ¹³C in Amazonian wood samples across 47 sites (circles), showing spatial variation in predicted values (‰) based on the final Random Forest model. (b) Map of associated predictive uncertainty expressed as the standard deviation (SD) of predicted δ¹³C values. Darker red areas indicate regions with higher model uncertainty.

3.2.2 δ¹⁵N isoscape

For δ¹⁵N, VSURF initially selected total soil nitrogen (r.nitrogen), soil organic carbon (r.soc), mean annual net primary productivity (r.npp), fire frequency (r.fire), atmospheric deposition of volcanic material (r.volc), and three geological age metrics (r.minage_geol, r.meanage_geol, and r.maxage_geol). After removing redundant predictors, only r.meanage_geol was retained among the geological metrics (r>0.98), and only r.soc was kept relative to r.nitrogen due to collinearity. Both r.fire and r.volc contributed negligible predictive power and were excluded to improve parsimony and ecological interpretability.

The final model therefore used mean geological age, soil organic carbon (SOC), and net primary productivity (NPP) as predictors. Among them, r.meanage_geol contributed most to node purity, followed by r.soc and r.npp. Cross-validation using 44 sites indicated that the Random Forest explained 67 % of the variance in δ¹⁵N, with RMSE = 1.10 ‰ and MAE = 0.93 ‰. δ¹⁵N measurements were unavailable for three sites due to insufficient nitrogen signal.

The scatter plot of observed versus predicted δ¹⁵N values demonstrates a predominantly linear trend across the entire range of variation (Fig. 7a). A systematic deviation is apparent, consistent with the conservative smoothing observed for δ¹³C, with higher δ¹⁵N values being underestimated and lower values being overestimated. The majority of data points are densely populated in the intermediate sector, specifically between 3.5 ‰ and 5.5 ‰, whereas fewer observations are located at the extremes of the distribution.

Spatially, the residuals are distributed heterogeneously across the Amazon, with both positive and negative deviations dispersed among the sampling sites (Fig. S5). As observed for δ¹³C, this heterogeneous spatial distribution implies the absence of localized bias, thereby supporting the robustness of the δ¹⁵N model on a regional scale. Partial dependence plots (Fig. 7b) reveal predominantly nonlinear responses. δ¹⁵N increases with mean geological age up to ∼1000 Ma, stabilizes between 1000–2000 Ma, and slightly decreases thereafter. With increasing SOC, δ¹⁵N rises up to ∼250 g kg⁻¹, declining at higher values. NPP shows a positive effect up to ∼2100 gC m⁻² yr⁻¹, followed by stabilization or a slight decrease.

The interaction between geological age and SOC (Fig. 7c) shows the highest predicted δ¹⁵N values (∼6 ‰) in very old substrates (>1500 Ma) with relatively low SOC (150–200 g kg⁻¹), and the lowest values (∼3 ‰) in younger substrates (<1000 Ma) with high SOC (>250 g kg⁻¹). This highlights the combined influence of geology and SOC as stronger than either predictor alone. The three predictor layers are provided in the Supplement for interpretation (Fig. S6).

Utilizing the final model calibrated with 44 sampled sites, a δ¹⁵N isoscape was generated for the Amazon Basin (Fig. 8a). Predicted values range from ∼2 ‰ to 9 ‰, with higher δ¹⁵N in the northern, southern, and central-eastern regions, and lower values predominantly in the southwestern, southern, and eastern Amazon, largely reflecting geological age gradients with additional SOC effects.

Predictive uncertainty (Fig. 8b) is higher in areas with sparse sampling, particularly in the northern and eastern basin. Substantial uncertainty also occurs in parts of the south despite relatively dense sampling, possibly reflecting local δ¹⁵N variability, limitations of the selected predictors, or unaccounted environmental processes influencing nitrogen dynamics.

Figure 7 Random Forest regression performance and environmental partial dependence plots (PDPs) for δ¹⁵N in wood. (a) Relationship between observed and predicted δ¹⁵N values, with a fitted regression line (black) and the 1:1 reference line (dotted). (b) PDPs of mean geological age of the bedrock (left), soil organic carbon content (center), and mean annual net primary productivity (right) on predicted δ¹⁵N values. Solid black lines represent the average predicted response when varying each predictor independently. Histograms below each plot indicate the distribution of observed values for the corresponding variable, providing context on sampling density. (c) Interaction surface between mean geological age and soil organic carbon content showing their combined effect on predicted δ¹⁵N. The color gradient represents predicted δ¹⁵N values (‰), with warmer colors indicating higher values and cooler colors indicating lower values. Black contour lines delineate zones of similar predicted δ¹⁵N values.

Figure 8 δ¹⁵N isoscape and associated uncertainty. (a) Predictive isoscape of δ¹⁵N in Amazonian wood samples across 44 sites (circles), showing spatial variation in predicted values (‰) based on the final Random Forest model. (b) Map of associated predictive uncertainty expressed as the standard deviation (SD) of predicted δ¹⁵N values. Darker red areas indicate regions with higher model uncertainty.

4 Discussion

This study provides the first δ¹³C and δ¹⁵N isoscapes for Amazonian wood, offering a spatially explicit framework to unravel the ecological and environmental factors shaping isotopic variability across the basin. Our results indicate that δ¹³C patterns are predominantly controlled by relative humidity and vapor pressure, reflecting climatic and physiological regulation of water-use efficiency, whereas δ¹⁵N distributions are primarily influenced by geological age and soil organic carbon, underscoring the strong edaphic control on nitrogen cycling. Together, these findings reveal coherent geographic structures in both isotopic systems and highlight their complementary potential for ecological interpretation and for developing robust, science-based tools for timber provenance across the Amazon.

4.1 Intra-site variability

δ¹³C and δ¹⁵N values in wood showed a large range at the Amazon-wide scale. The variation arises from the interplay between large-scale environmental controls, such as climate, water availability, soil composition, and locally operating ecological and physiological processes at the individual and species level (Leavitt and Roden, 2022; Savard and Siegwolf, 2022). Our dataset, which includes a broad taxonomic representation (Fig. S1), reflects part of the extraordinary tree diversity of the Amazon (ter Steege et al., 2023). This exceptional biodiversity, bringing together species with contrasting functional strategies, inherently amplifies isotopic variability at local scales. In this sense, the high intra-site variability we observed in both isotopes is not surprising and confirms the importance of intra- and interspecific controls on isotopic composition in Amazonian woods, as also emphasized by (Batista et al., 2025). For example, when comparing two sites with multiple genera collected (Feijó and Itapuã do Oeste; Fig. 4a), we found systematic intra-genus and inter-genus differences in carbon isotopes.

Intra-genus variability of δ¹³C was not negligible (Table S2). In Itapuã do Oeste, Handroanthus exhibited the lowest dispersion (0.46 ‰), whereas Cedrela and Dipteryx displayed 0.92 ‰ and 1.2 ‰, respectively. This indicates strong ecophysiological contributions to local variability (McCarroll and Loader, 2004; Cernusak et al., 2013; Cernusak and Ubierna, 2022), partly reflecting interspecific differences within genera, since individuals were identified only at the genus level. Such variance is likely explained by contrasting water-use efficiency (WUE) strategies (Ponton et al., 2001; Wittemann et al., 2024), driven by differences in stomatal conductance, ontogenetic stage, canopy position, and microhabitat conditions, which together influence photosynthesis and carbon isotope discrimination (Farquhar et al., 1982; Camargo and Marenco, 2011; Salmon et al., 2011; Van der Sleen et al., 2014; Brienen et al., 2017; Cernusak and Ubierna, 2022). Cedrela also shows anatomical plasticity in response to climatic variation (Ortega-Rodriguez et al., 2024), adjusting xylem traits under drier conditions (Rodríguez-Ramírez et al., 2022), with species-specific differences in resilience that help explain intra-genus variability.

From a functional perspective, the observed hierarchy (Handroanthus > Dipteryx > Cedrela) corresponds to differences in WUE and associated traits such as wood density and xylem anatomy (Chave et al., 2009; Fichtler and Worbes, 2012; Hu et al., 2024). Ecologically, Handroanthus and Dipteryx are typically emergent trees, whereas Cedrela generally occupies the upper canopy (Ribeiro et al., 1999). Taller trees in Central Amazonia exhibit greater stomatal density and tighter stomatal control to sustain water transport along longer hydraulic pathways (Camargo and Marenco, 2011; Brienen et al., 2017). These structural and ecological differences likely explain the consistently higher δ¹³C values and greater WUE in Handroanthus, reinforcing that structural and physiological traits jointly determine part of the isotopic composition of Amazonian wood (Batista et al., 2025).

For nitrogen isotopes, individuals from Itapuã do Oeste (5.66±1.19 ‰) consistently exhibited higher δ¹⁵N values than those from Feijó (2.25±1.67 ‰; Fig. 4b), highlighting strong inter-site differences (Nardoto et al., 2014). Intra-genus variability was generally lower in Itapuã do Oeste, except for Cedrela (SD ±1.66 ‰; Table S2), consistent with the high isotopic heterogeneity previously reported for terra firme forests (Nardoto et al., 2008). While δ¹⁵N values reveal clear regional patterns, local variability can be substantial, as individuals may exploit distinct nitrogen sources or respond to micro-edaphic and hydrological variation depending on age, canopy position, or successional group. Systematic inter-genus differences were also observed in Feijó, with Handroanthus presenting higher δ¹⁵N and Cedrela and Dipteryx lower values, likely reflecting edaphic conditions where organic N predominates and species explore contrasting nitrogen sources (Michelsen et al., 1996; Craine et al., 2015). Cedrela and Dipteryx appear to rely more on surface organic N, while Handroanthus accesses enriched mineral N pools. In Itapuã do Oeste, δ¹⁵N values were narrower, likely because the more “open” nitrogen cycle of highly weathered, well-drained soils promotes greater homogeneity in nitrogen source use (Martinelli et al., 1999). Notably, even Dipteryx (a nitrogen fixer) showed relatively high δ¹⁵N, consistent with evidence that most Amazonian legumes rely predominantly on soil nitrogen rather than atmospheric fixation (Högberg, 1997; Ometto et al., 2006; Nardoto et al., 2008, 2014).

Overall, δ¹³C variability was mainly structured at the taxonomic level, reflecting functional differences among genera, whereas δ¹⁵N variability was primarily driven by site-specific factors, with consistently higher values in Itapuã do Oeste than in Feijó. These patterns align with Martinelli et al. (2021), who analyzed foliar isotopic composition from evergreen forests across Brazilian biomes. Together, they indicate that high intra-site variability reduces the discriminatory power of carbon and nitrogen isotopes at the individual-tree scale, limiting their effectiveness for fine-scale provenancing.

At broader spatial scales, both Batista et al. (2025) and our study show that inter-site variance exceeds intra-site variability, reflecting the influence of large-scale environmental gradients. Building on these patterns, our results clarify the underlying drivers, showing that δ¹³C variability is primarily associated with climatic conditions, especially mean relative humidity, whereas δ¹⁵N reflects contrasting geological domains between Precambrian shields and Phanerozoic sedimentary basins (Fig. 3). This geographic structuring supports the potential use of δ¹³C and δ¹⁵N as complementary intrinsic tracers for applications in illegal wood provenance.

4.2 Relative humidity as the primary driver of δ¹³C in Amazonian wood

The δ¹³C dataset reveals a well-defined spatial gradient across the Amazon that follows patterns of relative humidity and WUE. Despite the ecological and climatic heterogeneity of the region, the Random Forest model performed robustly (R²≈0.60; RMSE ≈0.59 ‰) and identified mean annual relative humidity and mean annual vapor pressure as the main predictors (Fig. 5a, b). Isotopic trends broadly follow the moisture transport across the basin, with Atlantic vapor accumulating against the Andes (Salati et al., 1979). In the northwest, precipitation often exceeds 3000 mm and relative humidity approaches 90 %. Local vapor recycling accounts for up to 30 % of rainfall and is particularly intense in the west and southwest (Dominguez et al., 2022), reinforcing humid conditions (Shi et al., 2022) and favoring more negative δ¹³C values. In the south, southeast, and far north, precipitation is lower and relative humidity can drop to 40 %–60 % (Marengo et al., 2018; Espinoza et al., 2019; Marengo et al., 2024; Aprile et al., 2024), increasing dependence on evapotranspiration, climatic vulnerability, and resulting in less negative δ¹³C values.

In shield regions, precipitation is lower (1600–2200 mm), the dry season is prolonged, and relative humidity typically ranges between 40 %–60 % (Aprile et al., 2024), leading to less negative δ¹³C values. Roraima illustrates this contrast well: the drier north and east, dominated by savannas and open forests (Barbosa et al., 2007), show higher δ¹³C values, while the west and southwest, influenced by humid northwestern circulation, display more negative values (Barni et al., 2020). Within the Guiana Shield, Amapá also shows relatively less negative δ¹³C despite high annual rainfall, consistent with a more pronounced dry season and lower mean atmospheric saturation in the eastern Amazon (Marengo et al., 2001; Aprile et al., 2024), likely increasing vapor pressure deficit and aligning with our isotopic observations.

The relationship between δ¹³C and atmospheric vapor circulation reflects underlying ecophysiological processes. Vapor pressure deficit (VPD), defined as the difference between the saturation vapor pressure (e_max) and the actual vapor pressure of the air (e_a), is strongly correlated with relative humidity ($\mathrm{rH}=e_{\mathrm{a}}/e_{\max }$; $r\approx -\mathrm{0.96}$, Pearson), and is a key driver of stomatal conductance and δ¹³C variation in C₃ plants (Cernusak et al., 2013; Novick et al., 2016). In our model, mean annual relative humidity and vapor pressure together capture spatial gradients in atmospheric moisture, improving representation of δ¹³C discrimination through VPD-related processes. In humid regions such as the Solimões–Amazon corridor, floodplains and closed canopies maintain near-saturation conditions, reducing VPD, favoring stomatal opening, and enhancing discrimination against ¹³C, resulting in more negative δ¹³C values (Farquhar et al., 1989; Lloyd and Farquhar, 1994; Cernusak et al., 2013). In the south and southeast, higher temperatures increase atmospheric moisture capacity, relative humidity drops sharply in the dry season, evaporative demand intensifies, stomatal limitation increases, and δ¹³C values become less negative.

Partial dependence plots reflect this complexity: δ¹³C increases under intermediate conditions of relative humidity (∼63 %–68 %) and vapor pressure (∼2.7–2.85 kPa), when evaporative demand is intense, but decreases again under highly humid conditions (r.hurs_mean >70 % and r.vapor >2.85 kPa), when the atmosphere approaches saturation and discrimination increases (Fig. 5b, c). The interaction surface reinforces these patterns, with maxima under intermediate conditions and minima under saturated conditions, particularly along the humid Solimões–Amazon corridor. More negative values were also observed under low humidity and vapor conditions, which lack physiological support from the $c_i/c_a$ relationship (Farquhar et al., 1989) and may be associated with the lower density of observations in these ranges of the gradient (Fig. 5b).

Modeling uncertainty is highest in the east, north, and transition zones between the humid Solimões–Amazon corridor and less saturated areas (Fig. 6b), reflecting sharp humidity and vapor gradients (Fig. S4). Along these boundaries, small atmospheric changes may generate variable physiological responses, increasing predictive dispersion. Conversely, the southern and southeastern Amazon exhibit lower uncertainty, consistent with more homogeneous hydroclimatic conditions well captured by the predictors. Uncertainty is also influenced by sparse sampling in the east and north and by intra-site taxonomic variability, reinforcing the need to expand sampling, particularly in transition regions, with additional sites, more replicates, or more taxonomically focused designs.

4.3 Edaphic factors as the primary driver of δ¹⁵N in Amazonian wood

The spatial distribution of δ¹⁵N is not predominantly influenced by climatic gradients, but rather by edaphic factors. Variations between sedimentary and crystalline regions, associated with organic matter dynamics and primary productivity, are captured by three principal predictors: mean geological age, soil organic carbon, and net primary productivity. Together, these emerged as the strongest descriptors of the observed δ¹⁵N patterns, with the Random Forest model performing robustly (R²=0.67; RMSE = 1.10 ‰; Fig. 7). Similar controls were identified by Sena-Souza et al. (2020), who found SOC and NPP to be among the most significant predictors of continental-scale soil δ¹⁵N across South America, alongside edaphic and climatic variables, highlighting the critical role of these factors in shaping isotopic patterns associated with nitrogen cycling.

In our model, mean geological age was the most influential predictor (Fig. S6a), closely mirroring the δ¹⁵N isoscape (Fig. 8a). The Amazon is structured into two major geological domains: ancient crystalline shields and younger Cenozoic sedimentary regions (Martinelli et al., 2025). Partial dependence plots show a gradual increase in δ¹⁵N with geological age, stabilizing at higher values above ∼1500 Ma (Fig. 7b). This aligns with evidence that old cratonic surfaces host deeply weathered soils (Ferralsols and Acrisols), whereas younger sedimentary areas sustain less weathered Cambisols and Fluvisols (Quesada et al., 2011). Higher foliar δ¹⁵N in clay-rich, highly weathered Oxisols (Ferralsols), where iron- and aluminum-oxide-rich soils favor denitrification and fractionating nitrogen losses, has also been documented (Nardoto et al., 2008). Together, these lines of evidence support the view that ancient, nutrient-depleted soils consistently sustain higher δ¹⁵N than younger, less developed sedimentary lowlands. Remarkably, similar trajectories have only been described in short-range soil chronosequences from the Hawaiian Islands, where foliar and soil δ¹⁵N values increase with substrate age as nitrogen availability peaks and fractionating losses accumulate, before stabilizing in very old, phosphorus-limited systems (Crews et al., 1995; Vitousek et al., 1995; Vitousek and Farrington, 1997). Our findings extend this pattern to a regional scale, demonstrating for the first time a consistent association between geological age and δ¹⁵N across the Brazilian Amazon.

This Amazonian dichotomy likely reflects contrasts in bedrock porosity, mineralogy, and soil type that translate into differences in water-holding capacity, drainage, and soil redox conditions (Huscroft et al., 2018). Younger sedimentary lowlands, often influenced by Andean-derived or alluvial deposits, experience shallow water tables, frequent flooding, and pronounced hydrological fluctuations. These conditions favor denitrification and leaching but are offset by strong hydrological connectivity and continuous inputs of allochthonous nitrogen from upstream floods and sediments, which replenish the soil pool with isotopically lighter nitrogen sources, including substantial contributions from biological N₂ fixation (Martinelli et al., 1999; Hedin et al., 2009; Nardoto et al., 2008).

Conversely, ancient crystalline regions are well-drained, less prone to saturation, and receive minimal external nitrogen inputs (Gleeson et al., 2014; Huscroft et al., 2018). Over long pedogenic timescales, internally recycled nitrogen becomes progressively enriched in ¹⁵N through fractionating losses (NH₃ volatilization and denitrification), leading to consistently higher δ¹⁵N values (Martinelli et al., 1999; Davidson et al., 2007; Nardoto et al., 2008, 2014), a pattern clearly captured in our δ¹⁵N isoscape.

Amazonian soils are highly heterogeneous (Richter and Babbar, 1991; Quesada et al., 2011) and strongly shape biodiversity and ecosystem functioning (Cámara-Leret et al., 2017; Figueiredo et al., 2018; Schaefer et al., 2008; Tuomisto et al., 2003). Classical geochemical zonation distinguishes western Amazonia with young, nutrient-rich soils; central Amazonia with highly weathered, nutrient-poor soils; and northern/southern peripheries with intermediate cratonic-derived soils (Fittkau et al., 1975; Zuquim et al., 2023). Our δ¹⁵N isoscape broadly follows this pattern, but also reveals deviations. Central Amazonia showed relatively low δ¹⁵N values, closer to Andean-influenced regions than to cratonic areas, reinforcing the west–east δ¹⁵N gradient and the moderating role of hydrology (Nardoto et al., 2008, 2014; Sena-Souza et al., 2020).

Although δ¹⁵N values generally increase in the east, a lower-δ¹⁵N band follows the Amazon River corridor, reflecting strong hydrological connectivity. However, sites such as Belterra in the Santarém region retain high δ¹⁵N due to elevated terrain, well-drained clay-rich soils, and comparatively drier conditions that favor mineralization and fractionating nitrogen losses (Ometto et al., 2006; Nardoto et al., 2014).

Soil organic carbon was the second most important predictor. High SOC levels, typically associated with younger and hydrologically active environments, delay mineralization and enhance denitrification and leaching, maintaining lower δ¹⁵N values (Amundson et al., 2003; Houlton et al., 2006; Pérez et al., 2006). Conversely, older, well-drained soils exhibit lower SOC, enhanced mineralization, and higher δ¹⁵N (Nardoto et al., 2008; Craine et al., 2015). This relationship is clearly reflected in the interaction surface (Fig. 7c), with higher δ¹⁵N under older bedrock and low SOC and lower δ¹⁵N in younger, carbon-rich soils.

The least important variable was Net Primary Productivity (NPP). Although its overall contribution to the model was small, it significantly reduced prediction error. Previous studies in the Amazon indicate that NPP increases in more fertile soils, particularly those enriched in phosphorus and with higher foliar nitrogen (Aragão et al., 2009; Quesada et al., 2011). In such conditions, greater nutrient availability promotes vegetation growth and aboveground biomass accumulation, indirectly influencing nitrogen cycling. This effect was apparent only at a few sites, which explains the minor importance of NPP in the overall model.

The model shows a high degree of uncertainty in the Guiana Shield, reflecting the scarcity of samples and poor training in those regions. Elevated uncertainty persisted in Rondônia in the south despite good coverage, likely because the combination of contrasting soils and geomorphological heterogeneity generates hydrological microenvironments that increase isotopic variability. This pattern is consistent with the isoscape (Fig. 8a), which shows a mosaic of high and low δ¹⁵N values. In contrast, the south and southeast exhibit low uncertainty, consistent with more homogeneous environmental conditions. These results highlight the importance of expanding the sampling network, with priority given to critical regions such as the cratons, rock–sediment transition zones, and the Andean foreland. As with δ¹³C, this expansion should include a greater number of sites and replicates per locality, as well as the selection of taxonomically related genera or species to reduce non-climatic variability.

4.4 Limitations and perspectives

Although a systematic model bias is observed, characterized by a tendency to overpredict low isotope values and underpredict high values, the overall performance of the Random Forest models remains robust at the basin scale. This conservative smoothing behavior reflects a regression-to-the-mean effect inherent to ensemble tree-based approaches and is amplified when data density is highest in the intermediate range and sparser toward the extremes. Importantly, this bias is not spatially structured and therefore does not compromise the geographic consistency of the predictions; instead, it stabilizes basin-scale outputs by reducing overfitting.

In this context, while the models demonstrated good predictive performance, several important limitations remain, largely associated with the uneven sampling network across the basin. Uncertainty is particularly high in the cratons, in rock–sediment transition zones, and along the Andean foreland, where isotopic gradients are poorly represented by the available data. Inter-specific differences may also inflate the apparent spatial variability, reinforcing the value of constructing or calibrating isoscapes at least at the genus level for specific applications.

To increase the predictive power of these isoscapes, it is essential to expand the representativeness of the Amazon through new sampling efforts, especially in regions that remain underrepresented. It is also important to consider the taxonomic proximity among species in the sampling design, which may help reduce model uncertainties. Building a collaborative network, together with standardized sampling protocols, is crucial to ensure that all studies target the same radial position, allowing for more consistent comparisons and stronger discussions.

Even with these limitations, it is important to apply the models in real-world contexts and to assess their usefulness, costs, and potential integration into existing monitoring and enforcement systems. Overall, the combination of C and N isoscapes with geospatial predictions represents a promising and cost-effective pathway to increase the robustness of provenance assignments.

5 Conclusions

For the first time at the scale of the entire Amazon basin, we mapped δ¹³C and δ¹⁵N variations in wood. Our results show that these isotopic systems are informative for ecological processes and potential to play a complementary role in timber traceability. Carbon isotopes primarily reflect broad climatic gradients, whereas nitrogen isotopes capture more localized edaphic processes. Despite substantial intra-site variance, the combination of these isotopes has potential as a preliminary screening method that could identify potentially suspicious timber by batches rather than individual trees.

Although on their own they do not allow provenance at the level of individual logging sites, these isotopes establish a solid scientific foundation for regional-scale monitoring. This approach is particularly relevant in the Amazon given its vast territorial extent and the logistical complexity of enforcement. This study also advances the perspective of multi-isotope integration in which low-cost carbon, nitrogen, and oxygen analyses can serve as first-line screening followed by more precise and expensive systems such as sulfur, strontium, or genomic tools for definitive verification.

By demonstrating the feasibility of basin-wide wood isoscapes and applying scalable machine-learning approaches, this work contributes to tropical forest monitoring and to understanding forest functioning in relation to biogeochemical processes. These predictive frameworks can be applied in other tropical regions and offer a pathway toward global applications. The practical demands of forest governance make the isoscapes presented here a valuable complement to certification systems and sustainable management, and provide a strong foundational approach in efforts to secure a more reliable legal timber trade and to combat illegal logging in the Amazon.