Isoscape of Oxygen Stable Isotopes in Woods of the Amazon

1. Introduction

Stable isotopes have been extensively utilized in both ecological and forensic sciences [7,32,33]. The isotopic composition of elements such as carbon, hydrogen, nitrogen, oxygen, and sulfur varies based on the climatic, geographic, and ecological characteristics of the environment. These spatial variations, referred to as isoscapes, can be modeled using statistical and machine learning techniques such as kriging and Random Forest [6,15,18,22,24,37]. Isoscapes have proven valuable in tracing the geographic origin of materials ranging from timber and wildlife to human remains and narcotics [19,26,30,31].

These environmental signals are absorbed by plants and incorporated into organic tissues through well-established physiological fractionation processes [1,36]. In wood, cellulose is preferentially analyzed because it preserves the isotopic signal formed during biosynthesis and provides a more consistent representation of environmental conditions than bulk wood. Among stable isotopes, δ¹⁸O in tree-ring cellulose is particularly informative because it integrates both the isotopic composition of source water and leaf-level evaporative enrichment during transpiration. Consequently, δ¹⁸O values in cellulose provide a robust proxy for local hydroclimatic conditions and, ultimately, geographic origin [17].

Although the field of isotopic spatial modeling is well established, isoscapes remain scarce in tropical regions, particularly in the Brazilian Amazon, where applications are urgently needed. In particular, forensic approaches to timber tracking could benefit substantially from spatial isotopic tools. Despite existing regulatory frameworks, traceability challenges persist in parts of the Amazon timber supply chain, highlighting the need for complementary origin-verification methods based on chemical or genetic signatures. Recent studies have demonstrated the potential of isotopic models for this purpose [10,29,33,39], yet efforts remain fragmented and geographically limited.

In this study, we develop a regional δ¹⁸O isoscape for Amazonian tree-ring cellulose to examine the extent to which large-scale hydroclimatic gradients structure isotopic variability across the basin. We hypothesize that δ¹⁸O patterns in tropical wood are primarily governed by broad atmospheric and moisture-related controls rather than by local spatial dependence, reflecting the integration of source-water isotopic composition and evaporative enrichment processes at the leaf level. By evaluating how different statistical modeling frameworks represent these environmental controls, we seek to clarify the mechanisms shaping δ¹⁸O spatial structure in highly heterogeneous tropical forests. Through this approach, we aim to advance the conceptual and methodological foundation for constructing robust isotopic surfaces in the Amazon, where complex hydroclimatic dynamics and ecological variability challenge conventional spatial modeling assumptions.

2. Material and Methods

2.1. Isotopic Data Collection

To construct the δ¹⁸O-isoscape, we sampled 387 trees across 25 locations throughout the Amazon Basin (Figure 1). Species of commercial importance were initially targeted. However, due to the challenges of sampling the same species across the Amazon, expanded the sampling strategy to include mature trees regardless of their commercial significance.

Representatives from 58 genera distributed across 24 botanical families were collected. These forest species represent different ecological groups, including pioneers, secondary, and climax or later-successional, occurring in mixed or dense ombrophilous forests and reflecting the phytosociology diversity of the Amazon region.

Samples were collected from the base of trees harvested under authorized forest management plans within the Amazon region, in accordance with CONAMA Resolution Nº 406 of February 2, 2009[14]. The minimum cutting diameter (MCD) was ≥50 cm, and discs with a minimum thickness of 5 cm were obtained. Each wood disc was processed radially from the pith to the bark, and 2 mm thick laths were prepared, with their length varying according to the diameter of the collected disc.

The laths were proportionally segmented along the radial at (0% (pith), 25%, and 50% (heartwood), 75% (heartwood-sapwood transition), 100% (sapwood up to the bark)) (Figure 2). Each segment approximately 5 to 10 growth rings, enabling characterization of the radial variation in isotopic signal.

The tree-ring laths were chemically treated to obtain α-cellulose following the protocols established by the Isotopic Ecology Laboratory – CENA/USP [4], which are based on previously published and widely applied methodologies [23,28]. Isotopic analyses were partially conducted at the Stable Isotope Laboratory of the University of California, Davis (USA), and at the Stable Isotope Center at São Paulo State University (UNESP), Botucatu, Brazil. An intercalibration procedure was performed between laboratories to ensure analytical consistency and standardization of the results.

For oxygen isotope analyses, internationally recognized USGS reference materials (USGS54, USGS55, USGS56, USGS90, and USGS91) were used for calibration and quality control. The samples were analyzed for ¹⁸O/¹⁶O ratios using an isotope ratio mass spectrometer (IRMS) coupled with an automated pyrolysis device. The δ values were expressed as the relative difference between the isotopic ratio of the sample and that of the international standard (VSMOW for ¹⁸O) [35], as shown in Equation 1.

δ = ((R_sample / R_standard) - 1) × 1000

(1)

2.2. Environmental Data

Environmental variables such as Mean Annual Temperature (MAT), Mean Annual Precipitation (MAP), Potential Evapotranspiration (PET), and Vapor Pressure Deficit (VPD) were obtained from the WorldClim platform, version 2.1 (2020).

The data were extracted at a spatial resolution of 10 minutes (~340 km²). Each dataset contains 12 GeoTiff (.tif) files, one for each month of the year, numbered from January (1) to December (12). The extraction of these files was performed using QGIS software.

The Digital Elevation Model (DEM) was also obtained from the WorldClim platform, with elevation data derived from the Shuttle Radar Topography Mission (SRTM). These files were extracted at the same spatial resolution of 10 minutes, with an elevation precision of 10 meters. The rasters were acquired using QGIS software. Atmospheric Pressure (PA) was calculated according to Equation 2.

P(h) = P₀ × (1 - h / 44330)^5.255

(2)

where P(h) is the atmospheric pressure at different elevations, P₀ is the sea-level pressure (1013.25 hPa), and h is the altitude in meters from the DEM.

Relative Humidity (RH) was obtained using Equation 3. For the calculation of RH, water vapor pressure (kPa) raster files were first extracted from the WorldClim platform, with a spatial resolution of 10 minutes. The extraction of these rasters was performed using QGIS software.

RH = (C_a / es(T)) × 100

(3)

where RH is the relative humidity (%), Cₐ is the vapor pressure obtained from the raster, and es(T) is the saturation vapor pressure, calculated using Equation 4.

es(T) = 0.6108 × exp ((17.24 × T) / (237.3 + T))

(4)

where es(T) is the saturation vapor pressure (kPa), 0.6108 is a constant expressed in kPa, exp is the exponential function, 17.24 represents the linear dependence on temperature, 237.3 + T accounts for the increase in vapor pressure with temperature, and T is temperature in degrees Celsius (ºC).

The δ¹⁸O values of precipitation (δ¹⁸O_ppt) were determined by interpolating precipitation data from The Global Network of Isotopes in Precipitation (GNIP) database using IsoriX, an R package for constructing isoscapes using Generalized Linear Mixed Models (GLMMs) [16]. This process enables the creation of precipitation isoscapes using annual mean δ¹⁸O, generating raster data and extracting δ¹⁸O values from precipitation for each sampling point based on latitude and longitude.

The GLMM model in IsoriX accounts for spatial and temporal variations in isotopes and allows the response variable to follow a distribution other than normal, incorporating both fixed and random effects (Equation 5).

δ¹⁸O(S_i) = β₀ + β₁DEM + β₂T + μ(S_i)

(5)

where δ¹⁸O(S_i) is the isotopic value at spatial location S_i, β₀ is the intercept, β₁ and β₂ are coefficients for the covariates DEM and temperature, and μ(S_i) is the spatially autocorrelated error term, modeled using the Matérn function.

All environmental variables were extracted for each sampling location using the WorldClim database (Table 1). In addition, the isotopic composition of precipitation (δ¹⁸Oₚ) was interpolated for each site using the IsoriX package [16], based on the geographic coordinates and WorldClim environmental predictors. These interpolated δ¹⁸Oₚ values were later used as covariates in the modeling of δ¹⁸O in wood cellulose. The spatial distribution of modeled δ¹⁸Oₚ values reflect large-scale atmospheric processes, such as Rayleigh distillation and temperature-driven fractionation along air-mass trajectories [20,38].

These parameters (Table 1) were considered due to the global relationship between latitude, longitude, and δ¹⁸O_ppt, allowing for the prediction of local δ¹⁸O values [11,12,40,43]. Additionally, large-scale patterns in atmospheric vapor transport also influence δ¹⁸O in precipitation.

2.3. Construction of the Isotopic Model

Multiple Linear Regression Model

The isotopic model was developed using the caret package in R [25]. First, the predictor variables were standardized: the mean was subtracted from individual observations, and the resulting values were divided by the standard deviation, allowing for the observation of effect size. Subsequently, a maximal model (Equation 6) was constructed, incorporating potential predictor variables for δ¹⁸O in the Amazon region.

δ¹⁸O = β₀ + β₁VPD + β₂RH + β₃PET + β₄DEM + β₅PA + β₆MAT + β₇MAP + β₈δ¹⁸O_ppt

(6)

From this maximal model, all possible combinations of candidate variables were tested using the dredge function from the MuMIn package [3]. This process resulted in 256 unique models, from which the best was selected based on Information Theory (IT), considering those with ΔAICc < 4 [9,21]. As a result, five models were retained.

To evaluate model performance, Leave-One-Out Cross-Validation (LOOCV) was employed [8]. This method involves excluding one test data point, building the model with the remaining data, and calculating error metrics such as R², RMSE, and MAE. With model selection complete, the model equations were applied to variable values from the geospatial raster grids to construct the isoscape using the simplest (parsimonious) model.

We tested the regression model using both individual tree-level δ¹⁸O values and site-averaged δ¹⁸O values to evaluate its performance and applicability for timber provenance. Our main objective in making this comparison was to address a persistent question in timber forensics. In real-world applications, one typically has only a single wood sample to work with, rather than multiple samples from the same site that can be averaged to reduce variability. Therefore, while site-average models are valuable for identifying regional isotopic patterns, individual-level models are more realistic and directly applicable for forensic provenance, where predictions must be made from single samples. However, this approach introduces greater uncertainty due to natural intra-site variability, which in turn tends to reduce the model’s coefficient of determination. It is thus important to address the trade-off between model precision and practical applicability, balancing statistical performance with the realities of forensic casework.

Finally, we estimated Moran’s I statistic for the residuals of the regression model to assess spatial autocorrelation. The result (p = 0.73) indicated no significant spatial structure, justifying the exclusion of kriging from the modeling approach.

2.4. Random Forest

Random Forest (RF) is a machine learning (ML) technique applied in modeling[27,41]. We predicted δ¹⁸O of cellulose using RF and tested different covariates to select the most accurate model. Recursive feature elimination (RFE) was applied to select representative variables using the caret package. RFE constructs a model with all variables, ranks their importance, and iteratively eliminates the least important. Model performance was evaluated using repeated 5-fold cross-validation (10 repetitions) via the trainControl function. Metrics such as RMSE and R²were obtained from this process.

Like for the multiple regression model, we trained Random Forest models using individual tree-level δ¹⁸O values and site-averaged δ¹⁸O values to evaluate model performance and applicability to timber provenance. For provenance applications specifically, the model trained on site-mean values was selected for isoscape prediction and geographic assignment due to several critical methodological considerations. However, we also trained a Random Forest model using individual tree-level δ¹⁸O values in order to evaluate the full isotopic variability, including within-site variation. This allowed us to assess how the inclusion of biological noise arising from factors such as species differences, microclimatic conditions, and tree age impacts model performance.

The RF model was trained using caret package [25], with a grid search for tuning hyperparameters (Table 2). Variable importance was assessed, and predictions were compared to observed values. Environmental rasters were used to predict δ¹⁸O values across the Amazon Basin, and an isoscape was created and plotted using ggplot2 package [42], with the Legal Amazon boundaries serving as a reference. To better understand the structure of isotopic variability, we trained two Random Forest models: one using individual δ¹⁸O values from each tree and another using site-level means.

Table 2. Hyperparameter for variable selection in the random forest model.

Hyperparameter	Description	Unit
Longitude	Geographic coordinate (X)	Decimal degrees (º)
Latitude	Geographic coordinate (Y)	Decimal degrees (º)
VPD	Vapor pressure deficit	KPa
RH	Relative humidity	%
PET	Potential evapotranspiration	mm
DEM	Digital elevation model	Meters (m)
MAT	Annual average temperature	ºC
MAP	Annual average precipitation	mm
δ18Oppt	Predictive model of δ18O in precipitation	‰
δ18O	Isotopic ratio of 18O	‰

Table 2. Hyperparameter for variable selection in the random forest model.

Hyperparameter	Description	Unit
Longitude	Geographic coordinate (X)	Decimal degrees (º)
Latitude	Geographic coordinate (Y)	Decimal degrees (º)
VPD	Vapor pressure deficit	KPa
RH	Relative humidity	%
PET	Potential evapotranspiration	mm
DEM	Digital elevation model	Meters (m)
MAT	Annual average temperature	ºC
MAP	Annual average precipitation	mm
δ18Oppt	Predictive model of δ18O in precipitation	‰
δ18O	Isotopic ratio of 18O	‰

Table 3. Models selected by information theory, ΔAICc<4.

	Models	ΔAICc<4
1	${\delta }^{18}O={\beta }_{1}DEM+{\beta }_{2}Isorix+{\beta }_{3}MAP+{\beta }_{4}MAT+{\beta }_{5}PET+{\beta }_{6}RH+{\beta }_{7}VDP$	4456
2	${\delta }^{18}O={\beta }_{1}DEM+{\beta }_{2}Isorix+{\beta }_{3}MAP+{\beta }_{4}MAT+{\beta }_{5}RH+{\beta }_{6}VPD$	4458
3	${\delta }^{18}O={\beta }_{1}DEM+{\beta }_{2}Isorix+{\beta }_{3}MAP+{\beta }_{4}PET+{\beta }_{5}RH+{\beta }_{6}VPD$	4459

Table 3. Models selected by information theory, ΔAICc<4.

	Models	ΔAICc<4
1	${\delta }^{18}O={\beta }_{1}DEM+{\beta }_{2}Isorix+{\beta }_{3}MAP+{\beta }_{4}MAT+{\beta }_{5}PET+{\beta }_{6}RH+{\beta }_{7}VDP$	4456
2	${\delta }^{18}O={\beta }_{1}DEM+{\beta }_{2}Isorix+{\beta }_{3}MAP+{\beta }_{4}MAT+{\beta }_{5}RH+{\beta }_{6}VPD$	4458
3	${\delta }^{18}O={\beta }_{1}DEM+{\beta }_{2}Isorix+{\beta }_{3}MAP+{\beta }_{4}PET+{\beta }_{5}RH+{\beta }_{6}VPD$	4459

3. Results

3.1. Multiple Linear Regression

The model with the lowest AIC included the following climatic variables: mat, vpd, vap, rh, and pet. The coefficients produced by this model are shown in Table 1. The model built on the site-average data achieved an R² of 0.70, indicating that these predictors account for a significant portion of the variance in δ¹⁸O of the cellulose. The MAE was 0.56‰, and the RMSE was 0.68‰ (Table 4). Although not shown, the residuals followed a normal distribution, with no discernible pattern between the predicted values and the residuals. For the entire data set, there was a decrease in the R² to 0.42, and an increase of the MAE and RMSE to 0.92‰ and 1.10‰, respectively (Table 4).

The isoscape prediction grid spans longitudes between -74.00 and -43.50, and latitudes between -16.67 and 4.50, with a total of 12.392 points. The spatial grid has 183 cells on the x-axis (longitude) and 127 on the y-axis (latitude), with each cell measuring approximately 0.16 degrees, corresponding to spatial resolution of about 18 km² per cell. Model predictions based on site-average data ranged from 22.7‰ to 35.0‰, with a median value of 26.1‰ and interquartile range from 25.1‰ to 27.7‰. Predictions using individual tree δ¹⁸O values yielded a similar distribution, ranging from 23.8‰ to 33.2‰, with the same median (26.1‰) and interquartile range (25.1‰ to 27.7‰). These predicted values align closely with the observed δ¹⁸O measurements, which had a median of 25.8‰, interquartile rage (24.8‰ to 26.7‰) and ranged from 22.1‰ to 30.0‰.

A marked spatial gradient in δ¹⁸O values is evident across the Amazon region, with values decreasing from east to west, particularly along a southeast-to-northwest axis (Figure 3A). Although Figure 3 presents the isoscape generated from site-average δ¹⁸O values, a remarkably similar spatial pattern was observed in the isoscape derived from individual tree-level data (Figure S2), reinforcing the robustness of the observed isotopic gradient across the region.

Figure 3. Average isoscapes of δ¹⁸O in cellulose across the Amazon region, modeled using (A) multiple linear regression (MLR) and (B) random forest (RF), based on site-averaged δ¹⁸O values. Black lines indicate the boundaries of the Amazon region. Warmer colors represent higher δ¹⁸O values, while cooler colors indicate lower values.

Figure 3. Average isoscapes of δ¹⁸O in cellulose across the Amazon region, modeled using (A) multiple linear regression (MLR) and (B) random forest (RF), based on site-averaged δ¹⁸O values. Black lines indicate the boundaries of the Amazon region. Warmer colors represent higher δ¹⁸O values, while cooler colors indicate lower values.

Figure 4. Comparison of δ¹⁸O predictions in tree-ring cellulose across the Amazon region from the multiple linear regression (MLR) and random forest (RF) models. (A) Scatterplot of δ¹⁸O predictions from both models. The dashed red line represents the 1:1 relationship. (B) Spatial map highlighting areas where the difference between MLR and RF predictions falls within the ±1‰ interval (white), and areas with differences exceeding ±1‰ (gray).

Figure 4. Comparison of δ¹⁸O predictions in tree-ring cellulose across the Amazon region from the multiple linear regression (MLR) and random forest (RF) models. (A) Scatterplot of δ¹⁸O predictions from both models. The dashed red line represents the 1:1 relationship. (B) Spatial map highlighting areas where the difference between MLR and RF predictions falls within the ±1‰ interval (white), and areas with differences exceeding ±1‰ (gray).

The standard deviation of the predictions using site-average δ¹⁸O values ranges from 0.72‰ to 1.9‰. Descriptive statistics indicate an interquartile range of 0.75‰, to 0.83‰, with a median of 0.78‰ and mean of 0.83 ‰. The standard deviation of the predictions using individual trees δ¹⁸O values were higher, as expected. The median was 1.17‰, the interquartile range from 1.16‰ to 1.19‰, and the range from 1.15‰ to 2.05‰. Areas of low prediction uncertainty were distributed across the area, indicating relatively high confidence across much of the Amazon region (Figures S3 and S4). However, the areas with higher standard deviations tend to coincide with regions where sampling density is lower (Figure 1). Notably, many of these higher-uncertainty areas fall within the so-called "arc of deforestation", a zone of intense land-use changes and ongoing illegal logging activity.

3.2. Random Forest

The random forest model trained on the site-average data, yielding the lowest MAE of 0.64‰, and RMSE of 0.77‰ and a coefficient of determination (R²) of 0.67, indicating that the model explains approximately 67% of the variance in the data (Table 4). Thus, such quality metrics are comparable with those produced by the multiple linear regression model. These metrics deteriorate when the whole dataset is used to train the model. The MAE and RMSE increased to 0.90‰ and 1.12‰, with a R² of 0.44 (Table 4). Likewise, these metrics were also equivalent to those produced by the multiple linear regression model when using individual tree data.

Both the site-average and individual-tree isoscapes produced by the RF model showed similar spatial patterns, characterized by a southeast-to-northwest gradient in δ¹⁸O values. The site-average RF isoscape is presented in Figure 3B, while the individual-tree RF isoscape appears separately in the Supplementary Material (Figure S3). Despite these similarities, an important distinction emerged when comparing the RF and MLR model outputs: the range of predicted values was considerably narrower in the RF isoscape than in the MLR-based one (Figure 3A). Notably, the greatest spatial differences between the two models were found in the eastern Amazon, in the so called “arc of deforestation” where the models diverged most markedly (Figure 3B).

The standard deviation of the predictions using QRF and the site average dataset had slightly higher values than the standard deviations estimated using the MLR model. The median value was 0.87‰, the Q1 and Q3 were 0.70‰ and 1.10‰, respectively; and the range was from 0.2‰ to 2.0‰. However, the spatial distribution of standard deviation values across the Amazon was different than the one estimated by the MLR models. While the standard deviation isoscape by MLR produced a smooth surface with gradual transition (Figures S3 and S4), the standard deviation isoscape by QRF produced patchy surfaces due to the differences in model structure between the two methods (Figures S5 and S6). The standard deviation isoscape produced by the QRF model exhibits a patchy spatial pattern due to its data-driven estimation of local prediction uncertainty, while the smoother standard deviation surface from the MLR model reflects its assumption of homoscedasticity and linear relationships between predictors and response.

4. Discussion

This study represents the first effort to construct a δ¹⁸O isoscape from tree-ring cellulose across the Brazilian Amazon, aiming to contribute to the development of a multi-isotopic model to trace the origin of timber and combat illegal logging in the region. We compared the performance of two spatial interpolation approaches: multiple linear regression (MLR), and random forest (RF), a machine learning technique capable of capturing complex, nonlinear interactions between environmental variables. We also tested two data handling approaches: one using the δ¹⁸O values from individual trees, and the other using site-averaged values.

Our main findings show that the multiple linear regression (MLR) model performed slightly better than the Random Forest (RF) model, suggesting that the δ¹⁸O signal is primarily linear (Table 4). The comparison also highlights a central conceptual distinction in isotopic modeling for timber provenance. Models trained on site-averaged values capture regional or population-level isotopic patterns, whereas models based on individual samples explicitly account for within-site variability and are theoretically more appropriate for forensic provenance, where predictions must be made for single trees. Understanding this theoretical distinction is essential when designing models for enforcement or ecological applications. Additionally, it is crucial to evaluate how the use of individual versus site-averaged data affects the resulting δ¹⁸O isoscapes, as this may influence their spatial resolution and interpretability.

Both methods yielded isoscapes that reveal similar broad spatial patterns, with lower δ¹⁸O values in the western Amazon, higher values toward the east, and intermediate values in the central region. A latitudinal trend was also apparent, with elevated values at northern and southern extremes and lower values in the central Amazon. These patterns suggest that, despite differences in their modeling approaches, both methods can capture some of the large-scale environmental gradients influencing δ¹⁸O distribution.

The east-to-west decrease in δ¹⁸O values of cellulose mirrors the isotopic pattern observed in rainfall across the Amazon Basin [1,37]. This trend likely reflects Rayleigh-type fractionation, whereby heavy oxygen isotopes (¹⁸O) are progressively removed from atmospheric moisture during condensation and precipitation events along the prevailing Easterly moisture pathway [2,13]. As air masses travel inland from the Atlantic Ocean, successive rainfall events lead to the gradual depletion of ¹⁸O in water vapor, a signal ultimately captured in plant cellulose. Higher δ¹⁸O values in the east suggest shorter vapor transport distances and limited moisture recycling, whereas lower values in the west may indicate more extensive continental recycling and greater isotopic distillation. This east-west gradient underscores the dominant influence of Atlantic-derived moisture and its progressive isotopic modification during inland transport.

This example illustrates how synoptic-scale isotope effects on precipitation become imprinted in sessile organisms such as trees. The consistency between rainfall isotopes and cellulose δ¹⁸O reinforces the value of tree-ring isotopes as a powerful tool for deciphering complex climatic and hydrological patterns across large spatial scales.

Although the general spatial trends captured by the MLR and RF models were similar, the RF model predicted a noticeably narrower range of δ¹⁸O values compared to the MLR model. This pattern likely reflects the limited ability of Random Forest to extrapolate relationships between predictor and response variables beyond the calibration range, resulting in more conservative predictions within the bounds of the training data. In contrast, the MLR model can extend predictions linearly outside this range, potentially capturing broader gradients. Expanding the dataset through future sampling particularly in underrepresented regions will provide an opportunity to test the strength and stability of the linear relationships estimated by the MLR model and to evaluate whether these patterns persist as spatial coverage improves.

Most predicted standard deviations in the isoscapes were below 1‰; however, certain localized areas exhibited notably higher values, as shown in the standard deviation isoscapes generated by the MLR model (Figures S2 and S3). Many of these high-uncertainty regions overlap with the arc of deforestation. This convergence of sparse sampling, elevated model uncertainty, and ongoing environmental disturbance highlights a critical gap in current monitoring efforts. Future sampling strategies and model improvements should focus on reducing uncertainties in these vulnerable areas to improve the accuracy of provenance assessments and reinforce the application of isotopic tools in environmental enforcement.

Considerable variation in δ¹⁸O was observed among individual trees within the same plots, likely reflecting species-specific physiological traits and localized micro-environmental differences [5]. This highlights the importance of extensive within-site sampling to capture this variability and improve model resolution. While the current δ¹⁸O isoscapes are effective at distinguishing trees from isotopically distinct regions such as the far eastern and far western Amazon they are less reliable in isotopically homogeneous areas, particularly in the eastern Amazon, where logging pressure is most intense.

The predicted spatial range of δ¹⁸O values across the basin is approximately 3‰, while the model uncertainty, expressed as RMSE, is around 0.7‰. This suggests that regional differentiation is robust only where isotopic gradients exceed model uncertainty, for example, between western and eastern sites, but becomes ambiguous within low-contrast areas of central and eastern Amazonia. These results highlight both the potential and limitations of using δ¹⁸O isoscapes for timber provenance: while broad-scale distinctions are feasible, fine-scale geographic assignments remain uncertain. Expanding the spatial coverage of sampling and integrating complementary isotopic tracers (e.g., δ²H or ⁸⁷Sr/⁸⁶Sr) will be essential to enhance discriminative power in isotopically homogeneous regions where enforcement needs are greatest.

5. Conclusions

This study provides the first regional δ¹⁸O isoscapes for Amazonian wood and offers important methodological insights for isotope-based timber provenance in tropical forests. Although spatial autocorrelation was not statistically detected in our dataset, the coherent east–west gradients observed in the predictive maps indicate that isotopic structure likely exists but was not fully captured at the current sampling scale. Both multiple linear regression and random forest models generated broadly similar spatial patterns and comparable predictive errors; however, differences in predicted isotopic range were more evident in regions undergoing rapid environmental transformation, such as the arc of deforestation, highlighting that model selection may influence spatial interpretation and assignment sensitivity. By explicitly comparing site-averaged and individual-tree approaches, we demonstrate a fundamental trade-off between statistical robustness and forensic realism: while site means improve model fit and correlation strength, individual-level data more accurately reflect real-world enforcement scenarios, despite greater within-site variability. These findings underscore that isoscape development in highly heterogeneous tropical forests must balance generalizability with operational applicability. Future advances will depend on expanding spatial sampling coverage and integrating ecophysiologically relevant covariates, including traits such as stomatal conductance, successional status, and wood density, which influence plant water use and isotopic discrimination processes. In this sense, Amazonian timber isoscapes should no longer be regarded as static climatic surfaces, but as dynamic eco-physiological representations of forest structure, function, and environmental change.