The Development and Use of Isoscapes to Determine the Geographical Origin of Quercus spp. in the United States

Abstract

The stable isotope ratios of oxygen, hydrogen, carbon and sulfur from extracted wood of 87 samples of oaks from the United States were analysed. Relationships with climate variables and the stable isotope ratios of the 69 training dataset samples were investigated to a monthly resolution using long-term monthly mean climate data from NASA and the University of East Anglia’s Climate Research Unit, in conjunction with forecast data for hydrogen and oxygen isotope ratios in precipitation. These relationships were used to construct model isoscapes for oxygen, hydrogen, carbon and sulfur for US oak with the aim of using them to forecast isotopic patterns in areas that were not sampled and predict values in samples not used to construct the models. The leading predictors for isoscape generation were oxygen isotope ratios in January precipitation for oak oxygen isotope ratios, hydrogen isotope ratios in July precipitation for oak hydrogen isotope ratios, water vapour in April for carbon isotope ratios, and reflected shortwave radiation in March in combination with sulfate concentration in May for oak sulfur isotopes. The generated isoscapes can be used to show regions an unknown sample may have originated from with a resolution dependent on the rarity of the stable isotope signature within the United States. The models were assessed using the data of 18 samples of georeferenced oak. The assessment found that 100% of oxygen, 94% of hydrogen, 78% of carbon, and 94% of sulfur isotope ratios in the 18 test dataset samples fell within two standard deviations of the isoscape models. Using the results of the isoscapes in combination found that there were 4/18 test samples which did not fall within two standard deviations of the four models, this is largely attributed to the lower predictive power of the carbon isoscape model in conjunction with high local variability in carbon isotope ratios in both the test and training data. The method by which this geographic origin method has been developed will be useful to combat illegal logging and to validate legal supply chains for the purpose of good practice due diligence.

1. Introduction

Oak is one of the most traded timbers around the globe. Over 50 million cubic metres of oak have been exported from China every year since 2006 [1]. Though oaks are not typically considered to be threatened, the rampant trade in oak coupled with the ubiquity and high demand of oak furniture has caused significant deforestation in areas of the Amur Basin (China and Russia) [1] and the Carpathian Mountains. In China, this deforestation has led to the introduction of the National Forest Protection Program in 2016 to regulate logging in natural forests. As Chinese oak becomes less available, more manufacturers are turning to using US oak in their products as they can be relatively more confident of its sustainability and legality. Even so, as manufacturing sectors in Asia, where anti-illegal logging laws still have much room for improvement, continue to produce the majority of oak products, such as furniture and flooring, it is helpful to construct a tool to be able to verify the origin claims to ensure that the Lacey Act [2] of the United States, European Union Timber Regulation (EUTR, Regulation (EU) No 995/2010) [3], or Australian Illegal Logging Prohibition Act (ILPA) [4] are not violated. All aforementioned laws place emphasis on operators to carry out due diligence to ensure the declared origin and species of timber imports are correct and legal. This project focusses on creating a tool to aid enforcement to verify the geographical origin of oak. This analytical tool can be used to mitigate the risk of false origin declarations.

A secondary aim of this project was to build a physical collection of timber that is available for any researcher to access. The development of reference libraries of timber is the foremost limiting factor in the advancement of timber tracking tools. Many projects have focused on few participants obtaining reference samples [5,6,7] with no long-term ongoing sampling planned, or sampling plans have been broad, but instructions were not followed well by samplers such as in the ITTO project [8]. In this project, FSC and the US Forest Service co-ordinated the collection of reference samples using their considerable networks. The USA has been an ideal staging ground for future projects as collecting timber samples is largely non-controversial, there is an abundance of experts thanks to the US Forest Service and FSC, and the country is of relatively similar size to some of the countries that may need to be protected from illegal logging in future projects. If this project proves successful, there is an ambition to use FSC staff to collect samples of timber on a much wider and more regular basis for placement into the wood collection (xylarium) at the Royal Botanic Gardens, Kew.

Stable isotope analysis is a highly suitable tool to verify the origin of timber and has been demonstrated as an effective method to track the origin of timber [5,6,7]. It has an extensive history of use in legal cases in the food industry as demonstrated in Camin et al. [9] and is recommended by EUTR [3] and is therefore a valuable forensic tool. The success of the tool can be attributed to the fact that stable isotopes follow predictable patterns that are the direct product of geographically dependent processes such as the water cycle [10]. This means that relatively few samples are needed to map a country/region [11] giving it an advantage over other technologies such as population genetics and trace elements which often require greater sampling. The trade-off with the low sampling can be lower geographic resolution, though this can be remedied to an extent by modelling data [12,13,14] using related geographic variables and relevant statistics such as regression or kriging. This paper assesses the spatial variability of oxygen, hydrogen, carbon and sulfur isotope ratios of oak in the US, puts forward some models for isoscapes (isotope maps) based on regression and assesses their effectiveness at predicting stable isotope ratios in oaks.

2. Materials

Sixty-nine samples of various oak species (Quercus spp.) were collected from 13 US states comprising the training dataset. GPS co-ordinates of the samples were taken in the field using devices available to samplers including smartphones and Garmins. Notes were made about the diameter of the tree the samples were taken from along with the date of sampling and species of oak the collector identified in a sample pro forma (

Table 1. Species/variety/type of oak of the 69 training dataset samples as recorded by collectors.

Recorded Variety	Quantity Sampled
Chestnut oak	2
Quercus alba	55
Quercus falcata v. pagodaefolia	1
Quercus garryana	2
Quercus michauxii	1
Quercus rubra	2
Quercus rubra/garryana	1
Quercus suber	2
Red oak	1
White oak	2

and

Table 2. Quantity of samples collected in each state.

State.	Quantity of Samples
California	6
Kentucky	24
Maine	4
Maryland	3
Mississippi	3
New York	3
Ohio	1
Pennsylvania	6
Tennessee	4
Virginia	1
Washington	3
West Virginia	2
Wisconsin	9

); up to three leaf samples per tree were taken from 34 of the samples and placed in silica gel to use if future DNA extraction is deemed necessary and as herbarium vouchers. In total, 47 of the 69 samples were cross-sections cut from recently felled trees which were taken from the felled-end of the tree in some cases and at the crown end of the tree in others. Cross-sections were between 2 to 4 cm thick (Figure 1). Of the 69 samples, 22 were obtained from living trees using an increment bore. Samples varied in size with cores ranging from 10–18 cm by 0.5 cm and cross sections varying from 13–86 cm in width. These 69 samples comprised the training dataset for the isoscape models. An additional 18 samples were collected during the period of July 2017 to March 2019 from known locations in the USA in accordance with the same standards as the 69 training dataset samples (Figure 2). These 18 samples comprise the test dataset and were used to assess how effective the models were at predicting stable isotope ratios of oaks. The selection criterion for stable isotope analyses was that the sample must contain at least eight annual growth rings, this was to gain a robust mean value across the whole sample and mitigate the effects of atypical years of growth on the overall stable isotope ratios of the sample. For analysis, no delineation was made between heartwood or sapwood.

3. Methods

It is considered that isotopic fractionation from secondary metabolites in heartwood and the presence of functioning metabolic compounds such as proteins and starches in sapwood may have some effects on the overall isotopic composition of a sample. In the field of dendrochronology, this is typically accounted for by focusing only on a pure substance such as alpha cellulose. Agroisolab established a different measurement method [6] that is geared towards maximising the geographic information that can be extracted from a sample of timber and mitigating against the issue of exchangeable hydrogen rather than focusing on alpha cellulose. The reason for this is that methods such as those posed by Loader et al. [20] or Brendell et al. [21] involve the use of acids in the presence of water which may increase the risk of exchanging hydrogen with hydrogen in the cellulose. The method used in this project has minimal risk of exchangeable hydrogen, removes polar and non-polar contaminants and enables the measurement of sulfur isotope ratios in the same preparation which would not be available in an alpha cellulose preparation. Finally, as this method has been standard for Agroisolab since the mid-2000s, it should permit the generated data to be comparable with samples of timber historically analysed by Agroisolab. The trade-off of using this method is that measurements are not directly comparable to data from studies where alpha cellulose has been isolated for analysis and that there may be some influence of non-cellulose compounds on the isotope ratios that are measured. Nevertheless, the influence of non-cellulose compounds was internally assessed by comparing the hydrogen isotope data extracted using the Boner [6] method to non-exchangeable hydrogen isotope data derived from the same samples using a method described by Cheung [22].

Samples were initially dried at 103 °C before grinding/drilling. Cross sections were drilled using a 7 mm drill bit to the depth of 10–13 cm from the bark to the core. Increment bore samples were coarsely ground using a cutting mill (Retsch SM100–Haan, Germany). All materials were subsequently milled into a fine powder using a ball-mill (Retsch MM220–Haan, Germany). The powder was extracted in a soxhlet apparatus over 6 h with non-polar (Dichloromethane) and polar (methanol) solvents which were then dried in a laboratory-type drying cabinet for at least 1 h. Finally, the samples were stored in air-tight sample vials and weighed for analysis.

To avoid equilibration or humidity effects in the oxygen and hydrogen analysis, the weighed-in samples were equilibrated overnight in desiccators with a defined humidity of 10.6%. Afterwards the samples were vacuum dried for at least 2 h.

4. Measurement:

Three-point calibration was used to ensure the robustness of the measurements. Samples were compared at the beginning, end, and between each measurement run to a laboratory internal reference standard. Measurements are also reported relative to an internationally defined standard; for hydrogen and oxygen isotope ratio analysis, Vienna Standard Mean Ocean Water (VSMOW) is used. For carbon isotope ratio analysis, Pee Dee Belemnite (PDB) is used. For sulfur isotope ratio analysis Canyon Diablo Troilite (CDT) is used. Measurements are expressed in delta notation in accordance with the following Equation (1):

$$\delta =\frac{{\left(heavy isotope/light isotope\right)}_{Sample}-{\left(heavy isotope/light isotope\right)}_{Standard}}{{\left(heavy isotope/light isotope\right)}_{Standard}}\times 1000$$

(1)

Measurements are reported in permille (‰) and were made in accordance with processes outlined in Boner et al. [6].

¹⁸O/¹⁶O and D/H measurement: The high-temperature application uses HT-PyrOH with a covalently bonded silicon carbide tube (patented by Agroisolab GmbH, Julich, Germany) filled with glassy carbon chips and coal powder. The working temperature for pyrolysis is typically >1530 °C. To gain a higher precision the isotopes are measured in a master/slave configuration with two Isotope Ratio Mass Spectrometers (IRMS) (Isoprime, Elementar, Cheadle, UK). Each IRMS only measures one isotope ratio; D/H or ¹⁸O/¹⁶O. This configuration provides excellent stability because the magnetic field and accelerating voltage remain constant on each IRMS.

¹³C/¹²C measurement: Elemental Analysis (EA3000, Eurovector, Milano, Italy) in combination with IRMS (Nu Horizon, NU-Instruments—Wrexham, Wales). The working temperatures are 1021 °C for oxidation and 600 °C for reduction. Reduction is carried out in the presence of copper.

³⁴S/³²S measurement: Elemental Analyser (EA3000, Eurovector, Milano, Italy) with IRMS (Isoprime, Cheadle, UK). A one tube combustion (oxidation and reduction in one tube) is used to solve issues caused by SO₃. Combustion water is directly trapped with magnesium perchlorate at the end of the tube. The working temperature is 1000 °C.

5. Statistical Analysis

Latitude and longitude information was obtained from the sample pro formas (field records). Most collectors recorded this information in decimal degrees. Where collectors recorded latitude and longitude in degrees/minutes/seconds format, google maps was used to convert this information into decimal degrees. A CSV (Comma Separated Values) file of the data was created and uploaded into QGIS 2.12.3-Lyon (QGIS Development Team, 2012. QGIS Geographic Information System. Open Source Geospatial Foundation) and a shapefile with the co-ordinate reference system WGS 84 EPSG:4326 was saved. Inverse distance weighting (IDW) to the power 2 maps were created using SAGA GIS 2.1.2 (Departments for Physical Geography, Hamburg and Göttingen, Germany). Inverse distance weighting was selected to obtain a gross overview of geographic trends in the data. Distance weighting to the power 2 gives a greater weight to a mapped value over a given distance than a lower order function and is useful at representing spatial data without over exaggerating any given point’s influence. Though IDW does not give information about how predictable data are, it is a fair representation of the raw data of the samples on a map that has not been manipulated in any overt way.

SAGA GIS 2.1.2 was also used to create variograms of the isotope data for each parameter. Variograms can be used to gauge the spatial prediction quality of data. Variograms assess the variance of the data with respect over a nominal distance. The aim of using this statistical measure was to gauge whether the parameters followed a predictable geographic trend and whether there was enough variation in the data to allow for geographic discrimination. The distance over which the variogram assessed the data was 44.5 degrees with a lag distance of 0.5 degrees. In total, 59 of the 69 samples fit the distance range of variogram assessment. Ten samples were not included due to the lag distance size of 0.5 degrees. This was done so the variogram assessed the broader spatial predictability of the data across the lower 48 States rather than the local variance over shorter distances.

5.1. Linear Regressions—Model Building

Understanding the variables that influence or are related to particular isotope patterns is crucial especially in the field of geographic provenancing using stable isotope ratios. With current technology, it is not feasible to collect geotagged reference samples from every tree in existence, therefore there is always a level of prediction used in stable isotope provenancing. Understanding what influences isotope ratios should allow for statistical models to fill in the gaps where reference samples have not been taken with a greater level of accuracy than not having any data at all. This may eventually alleviate the burden of collecting reference samples in locations that are of political sensitivity to an extent provided there is enough data to generate a robust model. However, other researchers have created isoscapes using similar tools such as co-kriging [14,23]. These types of models are suggested to be advantageous as they can be multi-variate; they can incorporate data from multiple co-variates. Co-kriging models generate a variable confidence interval across a mapped region typically assuming that uncertainty increases with distance from a sample point. This is a useful function in many circumstances, however, the application put forth in this report requires a single value for uncertainty as the functions used to query the models look for values that fall within range defined by a set confidence interval.

5.2. Climatic Parameter Selection

Oxygen in cellulose originates from two primary sources, water and carbon dioxide [24]. Hydrogen in cellulose comes from water via photosynthesis. Isotope-based studies in oaks have shown this water to be groundwater [25] in some cases, though if other water sources are available to a tree it is reasonable that they can be used too. Global monthly grids for oxygen isotope ratios in precipitation were produced by Bowen and Revenaugh [12] to 20′ resolution. The oxygen grids were loaded into QGIS 2.12.3-Lyon along with the data shapefile. Using the point sampling tool (Borys Jurgiel, version 0.5.1) data were extracted from the monthly hydrogen and oxygen isotope grids at the same location as the reference samples. Linear regressions were performed in Microsoft Excel 2016 (Microsoft Corporation, Redmond, WA, USA) and fit was assessed by R² and visual inspection of the scatter plots.

Carbon in cellulose originates from CO₂ via photosynthesis. Stomatal conductance and water stress in plants are often cited as rate limiters which may cause carbon isotope fractionation within leaves [26,27,28]. Pollutants and vapour pressure have also been documented as variables that can influence stomatal conductance [26]. It should be noted that trees are complex plants and variables that affect carbon ratios in leaves may not perfectly mirror carbon ratios in timber cellulose as metabolic fractionation may occur in the production of celluloses and lignins [29,30,31] and physical fractionations may occur through moving sugars throughout the tree from photosynthesis to storage. Monthly global grids were downloaded from various sources for various years. The predictor variables were rainfall, precipitation, vapour pressure, water vapour, outgoing longwave radiation, tropospheric CO₂ concentration, tropospheric ozone, sulfur dioxide, tropospheric NO₂ and tropospheric CO. Two grids were selected that are similar “precipitation” and “rainfall”, these grids have different geographic resolution and cover different parts of the globe. Resolutions of the predictor variable grids ranged from 0.1 × 0.1 degrees to 2 × 2.5 degrees. To minimise the effects of extreme climates in certain years, long-term average grids were produced by summing the monthly grids (e.g., every grid for January 2006 to 2016) and dividing by the quantity of monthly grids used. Please refer to the supplementary materials for more details on the data sources and spatial resolution of the data (Tables S1 and S2). Using the point sampling tool (Borys Jurgiel, version 0.5.1) data were extracted from the monthly grids at the same location as the reference samples. Linear regressions were performed in Microsoft Excel 2016 (Microsoft Corporation, Redmond, WA, USA) and fit was assessed by R². Significance of each parameter in the multiple regression was assessed by the P value, where P values exceeded 0.05 the introduction of the parameter into the multiple regression was deemed not significant. If the introduction of a secondary or tertiary parameter improved the corrected R² of the multi-variate model beyond that of any single linear regression, this information was used to construct an isoscape incorporating information from those parameters using the equation to predict from multiple regressions Y = Constant + (B1 × (X1)) + (B2 × (X2)) + … BnXn. Where Y is the predicted outcome (in this case the forecast isotope ratio in oak), the constant is the intercept, B1 to Bn are the predictor variable regression coefficients, and X1 to Xn are the predictor variables (in this case the grids of climate or environmental variables).

Sulfur in timber has been cited to originate from several sources including soil, atmospheric deposition (emissions and volcanoes), and appears to be related to proximity to the sea [32]. Commonly, sulfur isotope ratios are measured to assess the source of sulfur in a material [33,34]. Factors that may influence sulfur ratios within timber may be related to a tree’s ability to uptake sulfur via transpiration and accessibility of sulfur in the soil via erosion or atmospheric deposition. Monthly global grids were downloaded from various sources for various years. The predictor variables were rainfall, precipitation, water vapour, cloud water content, plant canopy surface water, tropospheric nitrogen dioxide, tropospheric sulfate column mass density, and reflected shortwave radiation. To minimise the effects of extreme climates in certain years, long-term mean grids were produced by summing the monthly grids (e.g., every grid for January 2004 to 2018) and dividing by the quantity of monthly grids used. Using the point sampling tool (Borys Jurgiel, version 0.5.1) data were extracted from the monthly grids at the same location as the reference samples. Linear regressions were performed in Microsoft Excel 2016 (Microsoft Corporation, Redmond, WA, USA) and fit was assessed by R².

5.3. Model Assessment

The models were assessed in two ways. First, the proportion of training dataset samples (n = 69) that fell within the 95% confidence interval of all four isoscape models was checked along with how far away the nearest predicted region was in the event samples were not within this range. Secondly, the proportion of test dataset samples (n = 18) that fell within the 95% confidence intervals of the models individually, and in combination of the four models was assessed. Where samples did not fall within the confidence interval range of the models, the distance to the nearest region where the samples did match the confidence interval criteria was noted.

6. Results

6.1. Oxygen Isotope Ratios—¹⁸O/¹⁶O

The oxygen isotope ratios of the oak samples show a strong north/south trend with enriched oxygen isotope ratios of up to 26‰ in states that are geographically in the southern part of the USA such as California and Mississippi and depleted ratios as low as 23.2‰ in geographically northern states such as New York and Washington (Figure 3 and

Table 3. Stable isotope ratios of the 69 training dataset sample by State. Vienna Standard Mean Ocean Water (VSMOW); Pee Dee Belemnite (PDB); Canyon Diablo Troilite (CDT).

State	Quantity of Samples	18O/16O vs. VSMOW (‰)			D/H vs. VSMOW (‰)			13C/12C vs. PDB (‰)			34S/32S vs. CDT (‰)
		Mean	σ	Range	Mean	σ	Range	Mean	σ	Range	Mean	σ	Range
California	6	26.8	0.8	2.4	−98.5	6.6	14.0	−26.5	0.6	1.5	7.1	1.3	3.5
Kentucky	24	24.8	0.5	1.9	−85.7	6.6	25.6	−27.2	0.6	2.3	0.1	1.6	5.3
Maine	4	24.0	0.4	0.9	−96.9	1.7	3.7	−27.0	0.7	1.6	3.4	0.3	0.6
Maryland	3	24.2	0.5	1	−83.8	4.6	8.8	−26.3	1.2	2.4	0.7	0.3	0.5
Mississippi	3	26.4	0.7	1.2	−74.4	6.4	11.6	−28.8	0.7	1.4	2.4	0.7	1.2
New York	3	23.1	0.3	0.6	−101.1	1.3	2.5	−26.5	1.0	1.8	1.5	0.4	0.8
Ohio	1	23.8			−96.3			−26.5			−0.1
Pennsylvania	6	23.5	0.4	1	−93.7	2.0	5.9	−26.6	0.5	1.5	1.3	0.8	1.7
Tennessee	4	25.1	0.4	0.8	−81.0	3.1	7.6	−27.3	0.6	1.4	1.0	3.0	6.6
Virginia	1	23.4			−82.5			−27.3			0.9
Washington	3	21.9	0.5	0.9	−132.9	1.7	3.2	−25.7	0.4	0.8	4.0	0.3	0.5
West Virginia	2	25.1		0.2	−93.6		2.9	−26.0		1.5	0.2		1.3
Wisconsin	9	23.6	0.5	1.5	−92.0	1.9	5.9	−27.2	0.8	2.4	1.8	0.4	1.2

). The value range, 6.5‰, offers potentially a moderate degree of discrimination within the United States. The geographic trends are expected as precipitation isotope ratios show similar patterns to oxygen isotope ratios [12] and oaks assimilate water from precipitation into carbohydrates (e.g., sugars, cellulose). The trend is predictable as evidenced by the high variogram determination of 52.35% (

Table 4. Results of variogram analysis of the 69 oak samples across the United States. * Degrees of the Earth’s surface based on WGS 84 projection. At the equator, 1 degree is approximately 111 km.

Parameter	Variogram Determination (%)	Samples in Fitting Range	Lag Distance (Degrees *)	Fitting Range (Degrees *)
18O/16O	52.35	59	0.5	44.55
D/H	41.98	59	0.5	32.08
13C/12C	2.18	59	0.5	44.55
34S/32S	28.5	59	0.5	44.55

). As the oxygen in cellulose originates from CO₂ and H₂O, it is expected to show different geographic patterns to D/H ratios.

Overall, the best parameter for forecasting oxygen isotope ratios in US oak is the forecast oxygen isotope ratios in January precipitation [12] with an R² of 0.47 (Figure 4). Multiple regression assessment of January, February and December (or other combinations of months that were attempted) oxygen isotope precipitation grids to predict oxygen isotope composition in the 69 training dataset samples was not helpful or significant at improving the model as the P value of the secondary or tertiary predictors always exceeded 0.05. The 95% confidence interval for the oxygen isotope ratio model is +/−1.8‰, this value is high in relation to the work carried out by Gori et al. [14] where their 95% confidence interval ranged from +/−0.5‰ to +/−1.53‰ depending on the year and distance from the sample sites, however, it appears to be comparable to the maximum range of isotope ratios that may occur within a single oak tree ring which is given as 0.5–1.5‰ in Robertson et al. [35].

6.2. Hydrogen Isotope Ratios—D/H

The hydrogen isotope ratios of the oak samples show a different pattern to the oxygen ratios with elements of north/south and east/west trends (Figure 3B). The minimum value was −134.8‰ and the maximum value was −70.2‰ (range = 64.6‰) which should also afford a degree of in-country discrimination. The geographic pattern is predictable, with a variogram determination of 41.98%, this indicates that there are few unexpected ratios across the United States.

All R² values for monthly D/H isotope ratios in precipitation and D/H ratios in oak are quite high with the exception of January. This stands in contrast to oxygen isotope ratios where the highest R² value was in January. The highest R² value for hydrogen isotope ratios is in July (R² = 0.7, Figure 5). The 95% confidence interval for the created hydrogen isotope ratio model is +/−13.28‰. Again, this is higher than the values given for the models generated by Gori et al. [14] where their values for the 95% confidence interval ranged from +/−5.28‰ to +/−7.07‰ depending on the year and distance from the sample sites, but also appears to be consistent with the findings of the maximum circumferential variance of hydrogen isotope ratios in a single oak tree ring which is given as 5–20‰ in Robertson et al. [35].

6.3. Carbon Isotope Ratios—¹³C/¹²C

The carbon isotopes of the samples range from −25.3‰ in Northern regions such as Washington to around −29.5‰ in Mississippi (value range = 4.3‰). However, there is not such a clear north–south trend as demonstrated in Kentucky where the carbon ratios vary from −26.3 to −28.6‰, nearly the full range available across all the 69 reference samples. The narrow range in combination with the randomness of the values is well summarised by the low variogram determination of only 2.18%.

The best predictor of carbon isotope ratios of the climate/environmental variables that were assessed was water vapour pressure in April with an R² of 0.32 (Figure 6). The carbon isotope ratio model (Figure 7C) has a 95% confidence interval of +/−1.4‰. This value is relatively high, reflecting the uncertainty due to local variance seen in Kentucky. However, it is comparable to the maximum range of circumferential variance in a single tree ring which is given as 1–3‰ in Leavitt [11].

6.4. Sulfur Isotope Ratios—³⁴S/³²S (or δ³⁴S)

The sulfur isotope ratios of the oak samples follow a predictable pattern as indicated by the relatively high variogram determination of 28.5%. The lowest values are in Kentucky and Tennessee around −3‰ whereas the most enriched values are in California at +8.6‰. The range of values that exist across the reference samples (12‰) in conjunction with the fact there are not so many randomly different values across the United States suggests that sulfur isotope ratios should offer a degree of regional discrimination if more samples can be obtained.

Multiple regressions were performed to assess if the contribution of other parameters may improve the sulfur isotope ratio model. Combining Reflected Shortwave Radiation (RSWR) in March and precipitation in July a model with an R² of 0.69 was produced (adjusted R² 0.68). This is not a large improvement on simply using RSWR alone (R² = 0.66), though both parameters had P values below 0.05, therefore the multiple regression was significant. A further multiple regression was performed using RSWR in March and SO₄ concentration in March which produced an R² of 0.72 (adjusted R² 0.71).

6.5. Assessment of the Reliability of the Models

Of the 69 samples, 80% fall within the 95% confidence intervals of all four isoscapes. The remaining 14 samples show results between 18.56 km and 841 km away from the closest highlighted region. This error is demonstrated well in Figure 8C where the location of the sample is some distance away from the nearest region that fits with the 95% confidence interval of all the four models. Nevertheless, it should be stated that Figure 8C is able to show that regions in the same state (Washington) have signatures that are within the 95% confidence interval of the four models. Only four of the 14 samples have highlighted regions over 100 km away from their sampled location. In roughly half of the 14 outlier samples, the cause was sulfur isotope ratios that did not fit the sulfur isoscape model.

A total of 100% (18/18) of oxygen stable isotope measurements from the test dataset were within the 95% confidence interval range of the modelled oxygen isotope ratios at their location. Of the hydrogen stable isotope measurements from the test dataset, 94% (17/18) were within the 95% confidence interval range of the modelled hydrogen isotope ratios at their location. Of the carbon stable isotope measurements from the test dataset, 78% (14/18) were within the 95% confidence interval range of the modelled carbon isotope ratios at their location. Of sulfur stable isotope measurements from the test dataset, 94% (17/18) were within the 95% confidence interval range of the modelled sulfur isotope ratios at their location. Of test dataset samples, 78% (14/18) satisfied the 95% confidence interval range of all four isoscape models at their location. The 95% confidence interval range for the hydrogen, carbon and sulfur isotope ratios for sample 5 of 18 (taken near Roscommon, MI) were not satisfied, the nearest location where the criteria were satisfied was 600 km away in Wisconsin. This can mostly be attributed to a D/H isotope ratio that was more depleted than expected. If it were not for the depleted D/H isotope ratio measurement, the sample would have only been incorrectly assigned to be 54 km away. Samples 14, 15 and 18 of 18 had carbon stable isotope ratios that were outside of the 95% confidence interval of the carbon stable isotope ratio isoscape model. Samples 14 and 15 were both collected at the same site at the same site in near McCoole, MD. Four other samples of oak were collected at this site,

Table 5. Stable isotope ratios of the test dataset samples (n = 18) by State.

State	Quantity of Samples	18O/16O vs. VSMOW (‰)			D/H vs. VSMOW (‰)			13C/12C vs. PDB (‰)			34S/32S vs. CDT (‰)
		Mean	σ	Range	Mean	σ	Range	Mean	σ	Range	Mean	σ	Range
Maryland	6	24.3	0.2	0.6	−93.2	2.4	6.9	−26.2	1.3	3.5	0.2	1.0	2.6
Michigan	3	23.5	0.7	1.4	−103.5	9.6	18.4	−26.2	0.8	1.4	1.0	1.1	2.2
New Hampshire	2	23.0		0	−103.0		2.2	−26.9		0.7	2.0		1.3
Ohio	1	23.1			−90.7			−25.9			−2.4
Pennsylvania	2	24.4		0.5	−92.8		0.4	−26.4		1.1	2.3		0.9
Vermont	3	24.3	0.6	1.1	−100.6	7.9	13.8	−25.1	0.4	0.7	2.3	0.2	0.4
Virginia	1	24.8			−95.4			−27.3			3.1

shows that there was a 3.5‰ range in the ¹³C/¹²C stable isotope ratios at this site which is greater than the 2.8‰ range (95% CI) that the ¹³C/¹²C isoscape accounts for likely explaining the incorrectly predicted ¹³C/¹²C isotope ratio values for samples 14 and 15. Using the combination of models, the nearest locations where the 95% confidence intervals for samples 14 and 15 are satisfied are 45 km and 300 km away, respectively. Sample 18 of 18, collected near White River Junction, VT, had a ¹³C/¹²C isotope ratios that was 0.42‰ more enriched than expected at the location according to the ¹³C/¹²C isoscape. The nearest location where the stable isotope ratios of sample 18 fit within the 95% confidence intervals of the isoscape models approximately 100 km away near Middlebury, VT.

7. Discussion

7.1. Oxygen Isotope Ratios—¹⁸O/¹⁶O

Oxygen isotopic composition of precipitation in some months has a stronger relationship with the ¹⁸O/¹⁶O isotope ratios in the 69 training dataset samples than others. Figure 3 shows higher R² values during winter months than spring or summer months with the exception of July. On one hand, this is not as one may expect as photosynthesis and transpiration occur when the trees have leaves. On the other hand, autumn and winter generally produce more precipitation than summer in temperate forests. This could be due to the fact that precipitation across periods when the trees are not in leaf contribute towards the groundwater reservoir which the tree utilises in summer when in leaf. This may partially explain why the oxygen isotope ratios of the samples collected in California at the Institute of Forest Genetics are outliers; California drought may have led to a lack of groundwater and a need to irrigate the trees. As irrigation is unusual in silviculture, this may mean models would benefit if data from the Californian samples were removed (Figure 9). Nevertheless, the data for the Californian samples have been included in the models included in this paper. Furthermore, the influence of rainfall early in the season has been noted to break bud dormancy and trigger the growth of earlywood [36]. Therefore growth, and by relation stable isotope ratios, may be related to increases in soil water reserves.

July appears to be an important month in the summer for US oaks, no other summer month shows quite the R² as it. More information about the influence of each month may be garnered if isotope data from oaks around the world are added to the model. It should also be stated that no one month explained more than 50% of the variance in the oxygen isotope ratios in the 69 oaks; oxygen in cellulose/lignins originates from water and CO₂. Therefore, there may yet be further improvements that can be made to models.

7.2. Hydrogen Isotope Ratios—D/H

It is interesting that virtually all of the regressions of the D/H isotope data of the oak reference samples have relatively high values for R² with the precipitation isoscapes. It is also intriguing that the R² value for January is fairly low especially when contrasted with the fact that the R² value for January was one of the highest for oxygen isotope ratios in the oak reference samples. Perhaps this is a coincidence or perhaps it reflects that oxygen in cellulose originates from CO₂ and from H₂O whereas hydrogen in cellulose originates from H₂O. The highest single R² value for hydrogen isotope regressions was in July (R² = 0.7, Figure 5) which is a good basis to generate a model. Though the value for R² is high, 30% of the variance cannot be explained by the data from modelled precipitation isoscapes, meaning there is still room for improvement with the model. There could be several explanations for the model not being perfect ranging from model inaccuracy (e.g., the grid size of the Bowen isoscapes [12] is not fine enough to reflect the density of sampling or true variance of values in the US) to trees exchanging sugars via the mycelial network meaning the precipitation isoscapes are only part of the explanation.

7.3. Carbon Isotope Ratios—¹³C/¹²C

Based on the results of the regressions, water vapour in April, the month when oaks come back into leaf, is the leading predictive parameter for carbon ratios in US oaks based on the assessments made in this project (Figure 6), there may yet be further improvements to the model that can be made. Multiple regressions were carried out using the months with highest R² of other parameters. However, no additional variables caused a dramatic improvement to the model suggesting they are describing similar conditions (e.g., water vapour is related to vapour pressure and precipitation). Inclusion of pollution parameters (CO and NO₂) did not aid the predictive model either. A carbon isotope ratio model based on water vapour in April is available in Figure 7. This model predicts more extreme values than have been observed thus far with the reference data with areas within the great plains forecast to have significantly enriched carbon isotope ratios.

Pollutants such as Carbon Monoxide, CO₂, and NO₂, which have been stated to influence stomatal conductance [28,37,38] explain little of the variance in the carbon ratios of oak in this study. Tropospheric NO₂ in the US, particularly Eastern USA is relatively high compared to most regions on the planet; in July 2017 the AURA Satellite measured 108 billion molecules/mm² NO₂ in New York [39,40]. Only China has typically greater NO₂ concentration; in July 2017 the AURA satellite measured 178 billion molecules of NO₂/mm² [39]. Carbon monoxide has its greatest relationship with ¹³C/¹²C in oak in October (R² = 0.146, Figure 6), a time of year when oaks are typically shedding their last leaves. Most likely, this is a coincidence or the relationship refers to a covariate as there is no clear logical causal explanation for this phenomenon. The strongest relationship between NO₂ and carbon isotope ratios is in summer, particularly in July. In the context of the results from the oxygen and hydrogen isotope ratios, this may further suggest conditions in July are a key point in an oak’s calendar. To a degree, this appears to be consistent with the findings of Christison and Christison [41] that the growing months of leaf-shedding stress is confined to the months of June, July and August. This finding is also consistent with stomatal conductance experiments using pollutants [28], though it is noted that the relationship is not as strong as precipitation or vapour pressure which all have max R² relationships over 0.25. Various studies have emphasised the importance of pressure and vapour pressure’s influence on carbon isotope ratios [42]. Some have shown that carbon ratios may be related to radiation, however, the body of knowledge shows that any relationship is likely indirect as radiation can influence air pressure and water vapour pressure [26]. Water vapour pressure is related to stomatal conductance [26], it is also noted that low water vapour pressure can cause water stress in plants and can be a rate limiting step in the growth of trees and may also fractionate carbon isotope ratios [43].

Increased CO₂ emissions are expected to cause more depleted carbon ratios in trees [37]. Freyer and Belacy [38] demonstrated that across a 50-year span in tree rings, there is a 2‰ difference in carbon isotope ratios. The rationale for this is that as the carbon ratios of combustion carbon dioxide are more negative than natural CO₂, this is represented in tree rings. This study has not found this to be the case with the 69 oak reference samples. Instead, a very weak negative relationship between tropospheric CO₂ concentration and the carbon isotope ratios of oak was observed (Figure 6 and Figure 10). This may be due to the low resolution of the grid NASA has provided for CO₂ concentration (2⁰ × 2.5⁰ resolution) or due to the fact that the samples represent a ‘pooling’ of data from multiple tree rings in each sample. For perspective, the grid resolution is comparable to the area of the country of El Salvador. Though if this weak relationship is maintained when comparing data to higher resolution grids, this may be of benefit to timber tracking technologies as they may not need to be as concerned about the impact of increasing atmospheric carbon dioxide concentrations on the stable carbon isotopic composition of timber. Even so, perhaps this is not an ideal measure of this effect as several factors would provide a greater degree of information such as comparing the tropospheric carbon dioxide concentration in a single year with the corresponding growth ring in the tree as well as using carbon dioxide data that has a greater spatial resolution than 2⁰ × 2.5⁰.

7.4. Sulfur Isotope Ratios—³⁴S/³²S

To date, prevailing theories about sulfur isotope ratios in timber and agricultural products are that they may be related to plant or microbial fractionations [44], atmospheric deposition under certain conditions [45,46], or proximity to the sea [32,34,47]. Often, sulfur isotope ratio analysis is used just to identify the source of sulfur in geochemistry [33]. NO₂ concentration is related to pollution [48], March and November NO₂ concentrations appear to have some relationship with sulfur isotope ratios with R² of 0.3108 and 0.2724, respectively (Figure 11). This may suggest that atmospheric deposition has some influence on sulfur ratios, though the findings of Novak and Bottrell [46] must be noted; if sulfur has been in the soil for longer than sulfur from pollution (over 50 years) then the sulfur within the timber is predominantly natural. The seasonal patterns in R² and sulfur isotope ratios for cloud water content, precipitation, vapour pressure and water vapour may suggest that rates of transpiration play a role in the sulfur isotope ratio of oak. All the aforementioned parameters are related to transpiration [43] and the fact that strong relationships between the parameters and sulfur ratios exist while the trees are in leaf is not likely a coincidence.

The relationship of sulfur isotope ratios, sulfate and precipitation may be due to deposition of sulfate (SO₄²⁻) via precipitation; Mayer et al. [45] showed that the sulfur ratios of sulfate in precipitation are the same as sulfate in soil solution in aerobic forest soils. Using data from Reflected Shortwave Radiation and SO₄ together is currently the leading model for sulfur ratios of oak in the USA, the lack of a significant leap in improvement may suggest that both parameters are describing similar environmental attributes.

Consistent with the findings by Novak and Bottrell [46], it is doubted that the sulfur isotopes of the oak reference samples have been greatly influenced by pollution. The reason for this conclusion is that at the time of year when oaks are transpiring most (when they are in leaf), there are lower values of R² between tropospheric NO₂ concentration and sulfur ratios in the oak samples. Additionally, when values for NO₂ are included in multiple regressions with values for SO₄ or RSWR, NO₂ is not a significant parameter, its P value always exceeds 0.05. However, the impact of sulfate from precipitation may be an influencing factor as the addition of precipitation into a predictive model of sulfur ratios using multiple regression caused a small, but significant improvement in the model.

Finally, the sulfur isotope model appears to show similarity to the sulfur isoscape for human hair demonstrated in Valenzuela et al. [49]. In particular, Northern states such as Montana show similar patterns in sulfur ratios to the sulfur isoscape in Figure 3. However, there are some key differences between the patterns of the isoscapes in Ohio, Indiana, south Texas and Louisiana. This might be due to key differences between sulfur assimilation in humans and trees. As humans can consume food from anywhere in the world, Texans may have a different source of dietary sulfur than the sources of sulfur available to oaks.

8. Model Assessment Using the Test Dataset

The majority of test data samples (78% or 14/18) had isotope ratios that were in the 95% confidence interval of the model at their location. For the four samples that did not have values in the range forecast by the isoscape models, all were assigned to the USA but some were assigned up to 600 km away. This finding is important to consider as this error in assignment may be the difference in recommending further investigation is carried out on test samples. For sample 9, it appears that the resolution of the hydrogen isotope precipitation grid (approximately 18.5 × 18.5 km) that was used to predict the hydrogen isoscape was not able to account for the mountain location where sample 9 of 18 was collected. Perhaps including elevation, or some other high-resolution predictor grid as an additional predictor may help alleviate this problem in future. For samples 14, 15 and 18 of 18, the ¹³C/¹²C isoscape was unable to accurately predict their carbon isotope ratios within a 95% confidence interval. Again, this may be due to the resolution of the predictor grid (approximately 11 × 11 km) not being able to fully describe local variation in water vapour or the fact that ¹³C/¹²C isotopic composition of an oak is influenced by more variables than just water vapour. It also must be considered that the variogram assessment (Table 4) of the training dataset suggests that there is a weak spatial component to carbon isotope ratios. Table 3 and Table 5 also show that carbon isotope ratios have relatively broad ranges in sampling sites (e.g., 3.5‰ range in one site in Maryland n = 6) when carbon stable isotope measurements are taken across multiple growth rings. It is important to consider these findings in context of the natural range of ¹³C/¹²C isotope ratios in timber, which may only be approximately 10‰ in trees. Perhaps the limitation is more to do with the analytical method of measuring carbon isotope ratios across multiple growth rings, rather than carrying out individual ring measurements. Kagawa and Leavitt [50] produced models of carbon stable isotopes of Pinus edulis and Pinus monophylla based on a carbon isotope time series and individual ring measurements of ¹³C/¹²C that permitted the provenancing of these species with a precision of 114–304 km in southwestern United States. This method takes advantage of inter-ring variability in carbon whereas the method demonstrated here appears to be limited by it. Nevertheless, for those that do not have the ability to identify what growth ring in a piece of furniture or flooring corresponds to which year in a carbon isotope time series, measuring ¹³C/¹²C isotope ratios across multiple growth rings may be a more practical method despite its apparent limitations.

9. Conclusions

This project has produced detailed forecasts of oxygen, hydrogen, carbon and sulfur isotopes in the United States. The formulation of these forecasts lays the pathway for these techniques to be used to verify the origins of protected timbers. The production of these forecasts is potentially quite significant as it may provide a route to being able to verify the origin of timber in countries where it is not possible to obtain reference samples. It may also pave the way for forecasts of isotope ratios such as sulfur and carbon to be more routinely used in the verification of origin of other products such as fruits and vegetables. Over time, the ability to produce better and better forecasts will develop as the understanding of the causal factors grows.

The models demonstrated here are effective at predicting the stable isotope ratios of oaks from the USA, even considering that the test dataset samples were not collected in identical locations as training dataset samples. This finding is an encouraging demonstration of the predictive power of isoscapes. However, these models also have limitations. These include issues to do with the grid size of climatic predictors, the absence of including other predictor variables such precipitation amount-weighted precipitation isoscapes, elevation or others, and the highly locally variable nature of some parameters such as carbon isotope ratios. High local variability hampers large-scale prediction efforts. Perhaps other methods that make use of inter-annual variability in carbon isotope ratios between growth rings, such as the method proposed by Kagawa and Leavitt [50] may offer a better strategy than the carbon isotope ratio solution we put forth here if they can be used in practice on items such as furniture and flooring.

The ‘resolution’ of the analysis has been demonstrated to be variable and is based on the ‘rarity’ of a set of isotope ratios within a given geographic area. Nevertheless no strong inferences should be made about the resolution of isotope analysis using four stable isotope ratios from this paper as the resolution is shown to be variable dependent on: the rarity of the signature, the size of the grids used to make the forecasts, and the amount of available reference samples used to generate the model. It may even be possible to produce higher ‘resolution’ results with the same base data and more climatic information.

Finally, a map-based statistical method offers several advantages over pure-statistical models to verify the origin of timber. A map-based method may be able to integrate with other timber tracking tools to a better degree if they can be represented in map format. For example, it may be possible to use information from trace element analysis, chemical profiling, or Genomic analyses to further narrow-down areas that are forecast to be statistically consistent with the sample. Converting the statistical probability of a sample into a map gives both an ideal presentation tool and also may be a common ground for combining the results of various analytical methods.