Isotopic Characterization of Cedrela to Verify Species and Regional Provenance of Bolivian Timber

Illegal logging practices lead to loss of natural resources, community conflicts and economic problems. Illegal low-priced timber not only damages local and global forest markets (Abt Associates Inc. 2006; Blaser 2010), but also threatens sustainable management of tree species and forest areas (ABT 2017). The value of the wood trade associated with illegal logging in South American countries (Brazil, Colombia and Peru) has reached $682 million USD, with illegal logging rates ranging from 53 to 75% (Jianbang et al. 2016). The most common type of illegal timber trade is the fraud concerning false declarations of sites of origin based on falsified documents (Lowe et al. 2016). Increasing sustainability concerns call for developing tools to independently verify that the origins of timbers are in compliance with national and international regulations (Degen and Fladung 2007). A technique based on the intrinsic characteristics of the wood could help confirm the origin of the timber more reliably and thus help combat illegal timber trade.

The cellulose produced by trees to form wood tissue records the carbon and oxygen isotopic composition of the CO₂ and H₂O taken up by trees during photosynthesis and water uptake. These isotopic signatures can reveal important information about the environmental conditions the tree was experiencing at the time of growth. During photosynthesis, plants discriminate against ¹³C. The discrimination has been modelled by Farquhar et al. (1982) and is a function of the ratio of intercellular to ambient CO₂ concentration. This strongly correlates with water-use efficiency (Farquhar and Richards 1984) that can be estimated using tree-ring archives (van der Sleen et al. 2015). The oxygen isotopic composition (δ¹⁸O) in cellulose of tree rings records the isotopic composition of rainfall (Brienen et al. 2012). Recent studies have reported the relation between stable isotopes in Amazonian trees and precipitation (Evans and Schrag 2004; Jenkins 2009; Brienen et al. 2012; van der Sleen 2014). Stable isotopic composition of rainfall not only correlates to the amount of rain, it also shows distinct geographical patterns (Rozanski et al. 1993; West et al. 2010). Furthermore, both δ¹³C and δ¹⁸O are also correlated with elevation (Körner et al. 1988; Rozanski et al. 1993; Körner 1998), which could result in distinct isotope landscapes that can be useful to trace timber origin.

Stable isotopes have become an important tool to trace geographical origin of many biological materials. In particular, its application to ecology enables researchers to trace origin and migration of animals (Hobson 1999; Hobson and Wassenaar 1999; Rubenstein and Hobson 2004). Stable isotope analyses have also been successfully used in forensic applications. For example, stable isotopic composition of human hair is linked with geographic origin as it shows strong correlations with the geographic distribution of drinking water isotopic composition (Ehleringer et al. 2008).

In recent years, the analysis of stable isotopes has proven to be useful to accurately determine the origin of foods (Versini et al. 1997; Rossmann 2001; Oulhote et al. 2011; Driscoll et al. 2020) and drugs (Ehleringer et al. 2000; Booth et al. 2010). Geographic origin tracing by stable isotopes has also been tested for timber by Boner et al. (2007), Förstel et al. (2011), Förstel and Hützen (1983), Kagawa and Leavitt (2010), and Vlam et al. (2018). Although some of these studies have shown the potential for geographic identification at small scale (>14km) (Vlam et al. 2018; Watkinson et al. 2020), these studies reached different conclusions on the potential of timber tracing across larger spatial scales.

Stable isotopes have also shown promise for species identification (Chouvelon et al. 2014; Hamer et al. 2015), being successfully applied to fishes (Oliveira et al. 2011), birds (Militão et al. 2014), frogs (Dittrich et al. 2017) and archaeozoological material (Gorlova et al. 2015). To date, stable isotopes have not yet been used to identify tree species in a forensic context. Our study represents a first exploration on the potential of stable isotopes for the identification of Cedrela species and tracing back their site of origin.

Cedrela is an intensively harvested and traded tropical timber genus (Richter and Dallwitz 2000), widely distributed across climatic and environmental gradients throughout tropical America and the Caribbean islands (Wagenführ 2007), from moist to dry tropical forests and along a wide altitudinal range (Mostacedo et al. 2003; Navarro 2011; Navarro-Cerrillo et al. 2013). Cedrela has been traded as roundwood, lumber, veneer and exported predominantly to Europe and Northern America (Wagenführ 2007). Overexploitation of Cedrela species has caused population declines in recent decades (Mostacedo and Fredericksen 1999, 2001), motivating the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES) to include Cedrela species in its Appendix II (CITES 2019). Providing an objective, fast to implement, and reliable way to differentiate between Cedrela species and sites of origin across Bolivia, such as using isotopic compositions based on δ¹³C and δ¹⁸O, would greatly help in regulating and preserving these valuable and vulnerable tree species.

We expect isotopic signatures of carbon and oxygen to vary between sites, particularly at higher elevations, where spatial differences in terms of precipitation and topography are larger compared to the lowlands (Körner et al. 1988; Rozanski et al. 1993; Körner 1998; Brienen et al. 2012). Our research questions are: (1) Can we identify the geographic origin of Cedrela odorata and C. fissilis based on their isotopic signatures? (2) Can the relationship between precipitation, elevation and isotopic ratios be used to predict the likelihood of timber coming from specific sites? (3) Are the isotopic signatures different among Cedrela species? And, if so, to what extent do they differ? To address these questions, we collected Cedrela samples across multiple Bolivian forest types and assessed the isotopic variation within species in relation to precipitation, elevation and site of origin.

We performed species identification analysis on four Cedrela species (C. balansae, C. fissilis, C. odorata, and C. saltensis), and more-detailed origin identification on two of them (C. fissilis and C. odorata). Wood samples from 163 trees (2 wood cores per tree) were collected from eleven natural forest sites distributed across the species distribution areas in Bolivia (Figure 1). Cores were collected at breast height (130 cm) using a 5-mm-diameter increment borer (Haglöf Inc., Sweden). We targeted large/dominant trees in each site, as there is a height gradient associated with leaf δ¹³C values caused by the access to light and humidity (Ometto et al. 2002). Average tree diameter was 41 cm for both C. fissilis (range: 18–103 cm) and C. odorata (13–103 cm). Within each site, sampled trees were separated by at least 3 m for C. fissilis and 85 m for C. odorata. We tried to cover most of the species distributions in Bolivia and thus minimum distances between different sampling sites were 283 to 421 km for C. odorata, and 78 to 400 km for C. fissilis. We also sampled trees between sites, along roads and in villages, to improve spatial coverage and capture gradual spatial variation between sampling sites.

When species identification was not possible directly in the field, we collected botanical samples that were later identified by Alejandro Araujo Murakami, at the Museo de Historia Natural Noel Kempff Mercado (Bolivia). This was done for 27% of the sampled trees. Climate data were obtained from the reporting agency of Servicio Nacional de Meteorología e Hidrología (SENAMHI 2018; Table 1).

Samples were air dried and carefully polished beginning with ISO P40-grit (425–500 µm) sandpaper and using progressively finer sizes until ultimately sanding with ISO P600-grit (24.8–26.8 µm) (Orvis and Grissino-Mayer 2002). Cedrela species have clearly visible, ring-porous tree rings, with parenchyma bands as ring boundaries (Worbes 1999; Paredes-Villanueva et al. 2016). Tree-ring boundaries were identified using a binocular microscope (Leica S6E) coupled to an LED light source. Annual ring formation in Cedrela species is well established (Worbes 1999; Dünisch et al. 2002; Brienen and Zuidema 2005; Bräuning et al. 2009; Paredes-Villanueva et al. 2016). Tree rings were assigned to the calendar year in which their growth started, because growth season in these forests typically lasts from September of the current year to August of the following year (Schulman 1956). Next, annual ring widths within each increment core were measured, compared, and crossdated using TSAP/LINTAB software/hardware combination (Frank Rinn, Heidelberg, Germany) to a resolution of 0.01 mm. Missing and false rings, suggested by crossdated samples, were visually identified, re-measured and checked in an iterative process. We also estimated the mean first-order autocorrelation of raw tree-ring width data, a measure of the year-to-year growth similarity (AR1); mean intercorrelation among series, quantifying the similarity in width among trees (r); and mean sensitivity, which measures the year-to-year variability in width of consecutive rings (MS) for each site (Supplementary Material Table S1) using the dplR 1.7.2 (Bunn et al. 2021) R package. However, because the analyzed time-series spans were short (10 years), the results should be interpreted with caution.

We conducted exploratory analyses of δ¹³C and δ¹⁸O stable isotopes series of 10 years in bulk and annual samples separately. Wood flakes were sliced from the outermost 10 rings pooled into one sample (10-year bulk samples) of each core using a scalpel and microtome. Wood material from separate annual tree rings was also sampled from 3–4 trees per site (10% of the sample set) for all the sites, except Yapacaní. Crude cellulose was extracted following the adaptation of the Jayme-Wise method (Wieloch et al. 2011) described in Vlam et al. (2018) and then oven-dried at 50°C. For δ¹³C analysis, cellulose samples were combusted on an element analyzer, and for δ¹⁸O, cellulose was pyrolyzed using a glassy carbon reactor at 1450°C. Both were coupled to a continuous flow isotope ratio mass spectrometer (Sercon Hydra 20–20). The values of the isotopes ratios (δ¹³C and δ¹⁸O) per site were measured in parts per thousand (‰) according to the international Vienna Standard Mean Ocean Water (VSMOW) for oxygen, and Pee Dee Belemnite (PDB) standards for carbon.

Statistically significant differences between group means for species and sites were initially determined by one-way ANOVA. We then evaluated the relationship between isotopic composition, and local precipitation and elevation, using a mixed-effect modeling approach on annual measurements in C. fissilis, and on 10-year bulk measurements for C. fissilis and C. odorata with nlme 3.1–148 (Pinheiro et al. 2020) and MuMIn 1.43.17 (Bunn et al. 2021) R packages. We ran separate mixed models per species, to reflect potential species-specific responses. Sites were included in the models as random terms. In the case of annual time-series with multiple measurements per individual, tree identity was used as an additional random term. Annual precipitation of the last ten years was averaged for the studied sites. To be included in the 10-year bulk dataset, annual isotopic data were also averaged by the total years of measurement. We found no reliable climate data from Bajo Paraguá, and thus these samples were excluded from this analysis.

We performed a multi-step data analysis on the δ¹³C and δ¹⁸O data to evaluate whether this approach can consistently identify sites and species. As the dates of tree rings are generally not known in confiscated timber, we first needed to evaluate if taking into account annual variability may influence the main isotopic patterns and potentially hinder species or site identification. To test this, we compared two approaches for using the δ¹³C and δ¹⁸O annual data: (i) using raw averages of isotopic values, and (ii) using weighted averages based in their corresponding ring width. The idea behind the weighting correction is controlling for the fact that larger rings may contribute more to the mean isotopic value than small ones in the pooled sample. We performed a Pearson correlation to compare raw averaging and weighted averaging. After selecting the type of annual transformed data, these were then used as additional samples in the 10-year bulk dataset as input for the Kernel Discriminant Analysis (KDA). Pearson correlations between diameter at breast height (DBH) and stable isotope values were also calculated to assess if there was an ontogenetic effect in the sampled trees. Statistical analyses were performed in R version 3.4.3 (R Development Core Team 2017), using the ggplot2 2.2.1 (Wickham 2016) and ggpubr 0.1.6 (Kassambara 2017) packages.

The site differentiation potential was assessed by a discriminant analysis on the δ¹³C and δ¹⁸O in C. odorata (10-year bulk dataset) and C. fissilis (10-year bulk plus average annual data in one set). The Kernel Discriminant Analyses were done with the package ks 1.10.6 (Duong 2007, 2017) in R version 3.3.3 (R Development Core Team 2017), and were based on randomized samples and variables in every run. The classification results allowed us to assess the site differentiation measured by the frequency in which each site is assigned to either a correct or erroneous class.

KDA separates samples based on an a priori classification assignment (to specific site classes) and looks for the optimal non-linear combination of variables (component loadings) for maximal separation of the samples in a two dimensional space (Baudat and Anouar 2000). KDA's learning algorithm uses Bayes discriminant rule (Duong 2007) and needs to be trained in order to assess its discrimination power. Therefore, our data were split in two sets: 80% for training and 20% for testing the model in each of the analyses. Smoothed Cross Validation (SCV) error (Duong 2007) was applied to test the correctness of the site assignments. KDA testing generated confusion matrices showing the frequency at which each site was wrongly classified. After 100 randomization runs, we checked with which site a single sample origin could be confused most. The final classification error was expressed in percentage (%), where a cross-validation error of 0% would indicate that all samples were correctly assigned to their respective origins. Finally, the average errors per site identification across the 100 runs were obtained together with their corresponding standard deviations. For the species discrimination, KDA tests were performed on the 10-year bulk data by using four species as classes (C. balansae, C. fissilis, C. odorata and C. saltensis) and following the same steps described above.

Raw δ¹³C and δ¹⁸O data suggest differences between sites for C. fissilis and C. odorata (Figure 2). Higher isotopic values were observed in C. fissilis trees from Espejos for δ¹⁸O and Roboré for δ¹³C. C. odorata trees from Riberalta had the most contrasting isotopic signatures. Differences between other sites were smaller. Consistently, there were statistically significant differences between group means as determined by one-way ANOVA in δ¹⁸O (F_{5, 109} = 16.23, p < 0.001) and δ¹³C (F_{5, 109} = 8.17, p < 0.001) of C. fissilis sites; and in δ¹⁸O (F_{2, 27} = 5.82, p = 0.008) and δ¹³C (F_{2, 27} = 4.56, p = 0.02) of C. odorata.

The standard deviations for annual samples within trees ranged from 1.0 to 3.0‰ for δ¹⁸O (with averages between 22.0 and 25.3‰), and from 0.3 to 1.2‰ for δ¹³C (with averages between –27.3 and –26.0‰; Supplementary Material Table S2). Weighted and averaged annual isotopic data were significantly correlated for δ¹⁸O (Pearson correlation, r = 0.95, t₁₄ = 12, p < 0.001) and δ¹³C (Pearson correlation, r = 0.92, t₁₄ = 8.56, p < 0.001; Supplementary Material Figure S2). Moreover, both transformations followed similar isotopic patterns as the 10-year bulk samples (Supplementary Material Figure S3). Additionally, we found no significant correlation between tree-ring width and annual isotopic values for δ¹³C (Pearson correlation, r = –0.09, t₁₅₈ = –1.25, p = 0.21) nor for δ¹⁸O (Pearson correlation, r = –0.12, t₁₅₈ = –1.50, p = 0.13). Therefore, we decided to use the average annual values as additional samples within the 10-year bulk dataset for further analyses. The 10-year bulk isotopic values were clustered by site and species but had large overlaps (Figure 3).

The correlation between diameter at breast height (DBH) and stable isotopes (Supplementary Material Table S3) was not significant for C. odorata (δ¹⁸O Pearson correlation, r = –0.003, t₂₈ = –0.02, p = 0.99; δ¹³C Pearson correlation, r = 0.30, t₂₈ = 1.69, p = 0.10) or C. fissilis (δ¹⁸O Pearson correlation, r = –0.02, t₁₁₃ = –0.26, p = 0.79; δ¹³C Pearson correlation, r = –0.06, t₁₁₃ = –0.62, p = 0.53). This was also consistent within each site, with only Cobija, Riberalta and Roboré showing weak relationships between DBH and isotopes (Supplementary Material Table S3). Overall, any possible effect that size could have on isotopic ratios (Brienen et al. 2017) was hardly observed in these patterns.

Isotopic (δ¹³C and δ¹⁸O) composition in C. odorata showed no significant correlations with environmental conditions (mean site precipitation or elevation) in the 10-year bulk data. The relationship between stable isotopes in C. fissilis and mean precipitation and elevation varied depending on the spatial and temporal scale. At a larger temporal scale, i.e. when using 10-year bulk pooled data, we found no significant correlations between δ¹⁸O in C. fissilis and either precipitation or elevation, but δ¹³C had significant correlations with both variables (Table 2, Figure 4). Within-site linear model trends between elevation and isotopic (δ¹³C and δ¹⁸O) composition were not always significant and sometimes not consistent (Table 2, Supplementary Material Figures S1, S4, S5). It is important to note that elevation and precipitation values were highly correlated, and thus their significance when both variables were analyzed together in the model should be interpreted carefully. The variance inflation factor (VIF) of both variables in the models, however, were lower than the usually considered threshold of 5 (1.377 for 10-year bulk, C. fissilis and 1.424 for C. odorata), which suggested weak collinearity, and thus we decided to retain both variables and show the estimates of these models separately.

The classification matrix based on C. fissilis' 10-year bulk samples per tree (Table 3) estimated an overall site classification error of 60.2%. Yapacaní was the site with the least discriminating isotopic composition, with only 1% of samples correctly assigned to their corresponding origin and most of them wrongly assigned to Bajo Paraguá (82.5%). On the other hand, samples from Bajo Paraguá were mostly correctly identified (61.3%), and those misidentified were predominantly assigned to Espejos (11.8%). Roboré showed the highest classification accuracy (64.8%).

For C. odorata, mean total classification error was 70%. Riberalta samples were correctly assigned in 64.6% of cases, with some confusion with Cobija or Rurrenabaque (28.1–7.3%, Table 3). By contrast, the latter two sites showed correct assignments in less than 28% (Rurrenabaque) and 8% (Cobija) of samples.

There were statistically significant differences between species group means as determined by one-way ANOVA in δ¹⁸O (F_{3, 159} = 32.05, p < 0.001) and δ¹³C (F_{3, 159} = 28.57, p < 0.001). The KDA on the 10-year bulk samples showed clear species clustering, especially between C. fissilis and C. odorata (Figure 3a). Species discrimination for C. fissilis was remarkably accurate (95.3%, Table 4). Discriminating between the other three species was more challenging, with higher mean identification errors of 42.5% for C. odorata, 82.7% for C. saltensis, and 100% for C. balansae (Table 4A). Accuracy did not improve by only discriminating between the species with higher sample size (Table 4B).

In this study, we assessed isotopic differentiation in Cedrela species and their potential to trace wood geographical origin of C. odorata and C. fissilis. We hypothesized that δ¹³C and δ¹⁸O isotopic composition in Cedrela would reflect the spatial variation in environmental (McCarroll and Loader 2004) and topographic (Förstel and Hützen 1983) conditions, and thus could be used as a spatial proxy to identify timber origin. Our analyses show promising results but also highlight multiple challenges. Although stable isotopes showed significant correlations with precipitation and elevation at a large spatial scale (i.e. comparing average site conditions across sites), these relations were not always consistent within sites, suggesting that other environmental variables, rather than elevation or precipitation, also influence the variability in wood isotopic composition. This was also supported by our analyses of annual data, which showed no significant correlations between annual δ¹⁸O and δ¹³C values and precipitation (Table 2). Further analysis on site discrimination indicated a limited isotopic potential for tracing Cedrela timber geographic origin (Table 3).

Site classification accuracy was highly variable (between 1 and 64.8% for C. fissilis and 8.3 to 64.6% for C. odorata). This variation in classification accuracy may be explained by a high climatic overlap between some of the sites (similar precipitation conditions), which would also explain the misidentification between some of the sites. Alternatively, there could be species-specific signals or within-site variability in environmental conditions (e.g. in precipitation, elevation, light availability, or edaphic conditions) that may lead to the large differences in accuracy and overall low discrimination. Our discrimination analysis covered a wide spatial resolution, ranging from 3 m to 501 km for C. fissilis and from 85 m to 613 km for C. odorata samples. The overall low discrimination accuracy in some of the sites is consistent with what has been found for other tropical timber species based on δ¹³C and δ¹⁸O measurements, such as Erythrophleum in Cameroon and the Republic of Congo, where average accuracy was 35% (Vlam et al. 2018). Similarly, Vlam et al. (2018) reported a high variability in classification accuracy, from 46% up to 99%, at lower spatial scales ranging from 14 to 216 km. This suggests that the potential for implementing accurate and consistently reliable wood tracing based on these stable isotopes is likely site- and species-dependent. Also, our results highlight that the use of δ¹⁸O and δ¹³C alone may be insufficient to reliably account for variability in environmental conditions and scale.

Interannual variation in tree-ring δ¹⁸O and δ¹³C has been related to environmental conditions (McCarroll and Loader 2004; West et al. 2006). However, contrary to our expectations and to previous work (Dansgaard 1964; Gonfiantini et al. 2001), we found no significant correlations between annual isotopic data and precipitation or elevation (Table 2A). Low correlations between precipitation and δ¹³C and δ¹⁸O have also been reported in other tree species (Poussart and Schrag 2005; Cullen and Grierson 2007; Schollaen et al. 2013; Baker et al. 2015), and have been attributed to changing limiting factors of tree growth between years (Helle and Schleser 2004; Paredes-Villanueva et al. 2016). For example, because Cedrela is a relatively light-sensitive tree (Brienen and Zuidema 2006; Brienen et al. 2010), light availability might show a greater effect on tree development than precipitation in some sites (Brienen et al. 2010) but not in others. In theory, this periodical shift in the variable limiting growth could be taken into account by using annual isotopic data. However, in our study, using annual variation did not improve discrimination ability (mean error 70.1%), but rather added noise (Supplementary Material Table S4). Another source of uncertainty in our models may be the insufficient resolution and large scale of the environmental data (Kurita et al. 2009). Studying isotopes on a small scale while the limiting factor operates at a large scale can affect statistical analysis and consequent discrimination among sites. Indeed, large-scale environmental variation seemed to reflect site characteristics more reliably in our dataset when using 10-year pooled samples. The significant relationship between environmental conditions – such as precipitation and elevation – and the 10-year bulk isotopic signatures found on our exploratory analyses again reinforced the importance of considering the resolution and scale of environmental data. Analyzing longer time series of tree rings per tree could be a path forward to increase signals detection and accuracy of the results.

Although the potential of stable isotopes to trace timber geographical origin has previously been demonstrated in some species (Förstel and Hützen 1983; Boner et al. 2007; Kagawa and Leavitt 2010; Förstel et al. 2011; Vlam et al. 2018; Watkinson et al. 2020), timber tracing with a fine spatial resolution (i.e. higher precision) remains a challenge. We used a wide range of annual precipitation and elevation regimes to assess if their relation with trees will be mirrored in their isotopic composition (δ¹³C and δ¹⁸O). The overall inconsistencies between these environmental variables and our annual/10-year bulk isotopic data, in addition to the low discrimination among sites, suggest an influence by other governing factors at large spatiotemporal scales. For example, soil properties and El Niño-Southern Oscillation (ENSO), may not only influence isotopic composition in trees directly but also have marked effects depending on species and habitat. Previous research on δ¹³C and δ¹⁵N in coca (Erythroxylum coca) from South America found significant isotopic differences among sample sites and attributed it to differences in soil (for δ¹⁵N) and length of the wet season (for δ¹³C) (Ehleringer et al. 2000). Soil characteristics, particularly soil water availability and nutrients (Medina et al. 1995; Oliveira-Filho et al. 1998), have also been shown to have an impact on tree growth in dry forests (Toledo et al. 2011), and it is therefore likely to influence wood isotopic composition. This impact shows a spatial gradient (Murphy and Lugo 1986; Ceccon et al. 2006) that progressively shifts in wetter forests where light variation is more important for growth (Engelbrecht et al. 2007; Brienen et al. 2010). However, where soil chemicals are a limiting factor driving tree growth, it is important to further study which chemicals leave a signal in their wood composition.

The temporal-scale influence of ENSO is prevalent across South American forests, including on the isotopic composition of Cedrela species in Bolivia (Vuille and Werner 2005; Brienen et al. 2012; Baker et al. 2015). Vuille and Werner (2005) found higher δ¹⁸O values during the El Niño phase (drier) of ENSO, and lower δ¹⁸O values during La Niña events. Furthermore, at small spatial scale, precipitation was found to be associated with local circulation, and at a larger spatial scale, isotopic variations of precipitation were found to depend on the rain-out process in the surrounding region (Kurita et al. 2009). Although the rain-out history at a regional scale may be reflected in isotopic signal of precipitation, how stable isotopes composition varies in trees will depend on the site- and species-specific responses to these large-scale environmental factors (Helle and Schleser 2004; Poussart and Schrag 2005; van der Sleen et al. 2017). To our knowledge, temporal variability of stable isotopes has been studied for different species and sites in Bolivia (Brienen et al. 2012; Nijmeijer 2012; van der Sleen 2014; Baker et al. 2015; van der Sleen et al. 2015). However, for timber tracing, a good spatial representation is required. Although stable isotopes show promise, a better understanding of the between-species and within-sites variability is required before discrimination analysis can be widely and reliably implemented.

Our results showed distinct value ranges for each of the species' isotopic composition, with higher values for C. fissilis (ca. 24‰ δ¹⁸O, ca. –27‰ δ¹³C) compared to C. odorata (ca. 20‰δ¹⁸O, ca. –29‰ δ¹³C) (Figure 3). ¹³C discrimination is a function of the ratio of internal to ambient CO₂ concentrations (C_i:C_a) and water-use efficiency (Farquhar et al. 1982; Farquhar and Richards 1984). Under humid conditions C_i:C_a is larger resulting in more discrimination against ¹³C (more negative δ¹³C values). Under dry conditions discrimination is less, resulting in more positive δ¹³C values because of a relative enrichment in ¹³C. This was in line with our results, where the δ¹³C 10-year bulk values ranged from –31‰ for the moist-adapted species (C. odorata) to –24‰ for the dry-tolerant species (C. fissilis).

Although stable isotopes have been previously implemented to identify animal species (Oliveira et al. 2011; Chouvelon et al. 2014; Militão et al. 2014; Hamer et al. 2015; Dittrich et al. 2017), our study represents a first exploration on their potential use for taxonomic and forensics purposes in timber species. In our case, using 10-year bulk samples seemed to render highly accurate identification for C. fissilis with up to 95.3% accuracy, but seemed insufficient to reliably identify other Cedrela species, particularly C. balansae, which showed 100% error. However, the worst performing species in our classification analyses were also the least replicated ones, indicating the need for larger sample sizes for each species in future work. Consistent with our results, previous discrimination analyses using Direct Analysis in Real Time (DART) coupled to Time-of-Flight Mass Spectrometry (TOFMS) also found C. balansae challenging to identify, with the highest error of 46.1% compared to C. fissilis which presented the highest accuracy (91.3%) (Paredes-Villanueva et al. 2018). The current data, though limited, suggest that identifying C. balansae via isotopic signatures may require their use in combination with other (e.g. genetics, anatomical, Near Infrared Spectroscopy) approaches.

An additional source of identification error that should be taken into account concerns updates on the species classifications based on herbarium collections. Rapidly advancing tropical botanical studies mean that the morphological characteristics to tease apart Cedrela species are constantly being updated, which makes accurate identification via botanical samples sometimes challenging (Zapater et al. 2004; Koecke et al. 2015; Palacios et al. 2019). Species classification based on both morphological and genetic properties might help us to perform chemical discrimination for then well-established species. In addition, it will provide insights into the underlying causes of identification errors. These protocols can consequently be fine-tuned to discriminate between problematic species.

How much isotopic variation is there among all traded timber species? Are there any other species having the same isotopic values as any specific Cedrela species, which could increase the error rate for identification? Further work is needed to map the expected isotopic range for more timber species, in order to contextualize our results and provide widely applicable tools for sustainable timber trading.

Some considerations have to be taken into account to interpret our ¹³C and ¹⁸O results. First, a light-demanding species may differ in its isotopic composition if samples are collected at disturbed areas or if wood tissue was formed during clearing events (Brienen et al. 2010; van der Sleen et al. 2014). In our case, samples were collected in areas where there was previous intervention or use, as well as in dense forests with closed canopy. However, tree sensitivity and its isotopic composition may also be influenced by temporally shifting limiting factors (Helle and Schleser 2004) at a specific period. In case there is a possibility that the intervention in a forest or light availability plays an important role in the tree isotopic composition, we recommend taking into account the intervention degree in future discrimination analyses.

Second, we expected isotopic values to be reliable proxies for climate data (e.g. precipitation). Indeed, Cedrela trees' shallow roots would suggest a strong correlation between δ¹⁸O in trees and precipitation (Brienen et al. 2012). However, in our study, precipitation did not show high variation among C. fissilis sites (Table 1) nor significant correlation with the δ¹⁸O (Table 2, Figure 4). This suggests that precipitation may not always be the main driving factor of δ¹⁸O composition in Cedrela wood. Other possible oxygen sources that can have an effect on isotopic variation in plants include humidity, soil moisture, groundwater (Aggarwal et al. 2004), and clouds (Anchukaitis et al. 2008; Anchukaitis and Evans 2010).

Third, isotopic variation also depends on tree age and size (ontogenetic development phase) (Broadmeadow et al. 1992; McCarroll and Loader 2004; van der Sleen et al. 2015; Brienen et al. 2016). However, we found no supporting evidence for this effect influencing the discrimination potential of our isotope samples. Tree size (DBH) did not have a significant effect in the isotopic composition of our trees. We also found no significant correlation between DBH and stable isotopes in C. odorata, nor in C. fissilis. Similarly, removing trees with exceptionally small or large DBH (i.e. 2 C. odorata samples with <15 and >70 cm, and 4 C. fissilis samples with <20 and >70 cm) to minimize this potential size/age effects, had little effect on mean discrimination accuracy (+1.2% in C. fissilis and –0.5% in C. odorata).

In addition, during an analysis of timber with unknown origin, estimating the possible period of tree harvesting will be a minimum requirement to correct for the overall decrease in δ¹³C in tree biomass since the onset of the Industrial Revolution period (McCarroll and Loader 2004; van der Sleen et al. 2017). Finally, the low spatial variation in δ¹³C and δ¹⁸O could be improved by assessing multiple tree-isotope proxies or combining them with rare earth and trace elements (English et al. 2001; Kelly et al. 2005; Joebstl et al. 2010) to increase discrimination among sites.

Cedrela species are among the most important traded timbers species and have suffered from over-exploitation. We aimed at assessing whether timber isotopic composition (based on δ¹³C and δ¹⁸O) can be used to differentiate Cedrela species and sites of origin, providing a key regulatory tool for tropical forests. This approach could potentially improve the verification of timber origin certificates in accordance with regulations. However, our analyses reveal limited potential for discrimination of geographical origin using wood δ¹³C and δ¹⁸O. Sites assignment accuracy was low, at 30% for C. odorata and 39.8% for C. fissilis. This low geographical discrimination potential, using carbon and oxygen isotopic imprints alone, may result from multiple factors, including high climatic overlap between certain sites, high within-site variability, and lack of high-resolution environmental data. In addition, using precipitation and elevation as proxies to infer site of origin at a fine resolution (specific tree location) remains a challenge, given the high variability of isotopic composition (δ¹³C and δ¹⁸O) within trees. Yet, our results suggest that it is possible to distinguish between species based on stable isotopes. In particular, C. fissilis can be identified with a remarkably high accuracy (95.3%). We conclude that although much work is still needed to be able to use stable isotopes as a sole tool to identify tropical trees provenances, this methodology already shows potential for discrimination of some important timber species. Further developments in this direction will contribute to improve timber origin classification and thereby may become a tool in the international fight against illegal trade of (precious) tropical timbers.

We thank the Museo de Historia Natural Noel Kempff Mercado for their support on the species identification. This research was financed by the Dutch Organization for Internationalization in Education (NUFFIC) [NFP-PhD.14/61]. Fieldwork was logistically supported by Universidad Autónoma Gabriel René Moreno and financed by Alberta Mennega Stichting and The Rufford Foundation [18670-1]. This research was also supported by the International Foundation for Science (IFS), Stockholm, Sweden, through a grant to Kathelyn Paredes Villanueva [D/5914-1] for a second campaign of fieldwork and laboratory analyses. RDM was supported by the Swiss National Science Foundation through a Postdoc. Mobility Fellowship. We also thank the reviewers and editor of the journal who have provided valuable comments improving the quality of our manuscript.

Supplementary Material is available at http://www.treeringsociety.org/TRBTRR/TRBTRR.htm