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

: 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 reﬂected 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.


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 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. Forests 2020, 11, x FOR PEER REVIEW 3 of 23 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.   The green area represents a combination of the species distribution maps for Quercus spp. created by Little [15][16][17][18] and digitized by Thompson et al., [19].

Figure 2.
Locations of the training and test dataset samples in the United States. The green area represents a combination of the species distribution maps for Quercus spp. created by Little [15][16][17][18] and digitized by Thompson et al. [19].

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.

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): Measurements are reported in permille (% ) and were made in accordance with processes outlined in Boner et al. [6]. 18

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.

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.

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 2 and visual inspection of the scatter plots.
Carbon in cellulose originates from CO 2 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 2 concentration, tropospheric ozone, sulfur dioxide, tropospheric NO 2 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 2 . 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 2 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 2 .

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.

Oxygen Isotope Ratios-18 O/ 16 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). 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). As the oxygen in cellulose originates from CO 2 and H 2 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 2 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]. 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 2 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].    and forecast oxygen isotope ratios in monthly precipitation [12].

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 2 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 2 value was in January. The highest R 2 value for hydrogen isotope ratios is in July (R 2 = 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].

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 2 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 2 value was in January. The highest R 2 value for hydrogen isotope ratios is in July (R 2 = 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].

Carbon Isotope Ratios-13 C/ 12 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 2 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].

Sulfur Isotope Ratios-34 S/ 32 S (or δ 34 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 2 of 0.69 was produced (adjusted R 2 0.68). This is not a large improvement on simply using RSWR alone (R 2 = 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 SO4 concentration in March which produced an R 2 of 0.72 (adjusted R 2 0.71).

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.  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 2 of 0.69 was produced (adjusted R 2 0.68). This is not a large improvement on simply using RSWR alone (R 2 = 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 4 concentration in March which produced an R 2 of 0.72 (adjusted R 2 0.71).

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 shows that there was a 3.5‰ range in the 13 C/ 12 C stable isotope ratios at this site which is greater than the 2.8‰ range (95% CI) that the 13 C/ 12 C isoscape accounts for likely explaining the incorrectly predicted 13 C/ 12 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 13 C/ 12 C isotope ratios that was 0.42‰ more enriched than expected at the location according to the 13 C/ 12 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. 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 shows that there was a 3.5% range in the 13 C/ 12 C stable isotope ratios at this site which is greater than the 2.8% range (95% CI) that the 13 C/ 12 C isoscape accounts for likely explaining the incorrectly predicted 13 C/ 12 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 13 C/ 12 C isotope ratios that was 0.42% more enriched than expected at the location according to the 13 C/ 12 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.  Figure 3 shows higher R 2 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.

Carbon Isotope Ratios-13 C/ 12 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 2 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 NO2) 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 July appears to be an important month in the summer for US oaks, no other summer month shows quite the R 2 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 2 . Therefore, there may yet be further improvements that can be made to models.

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 2 with the precipitation isoscapes. It is also intriguing that the R 2 value for January is fairly low especially when contrasted with the fact that the R 2 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 2 and from H 2 O whereas hydrogen in cellulose originates from H 2 O. The highest single R 2 value for hydrogen isotope regressions was in July (R 2 = 0.7, Figure 5) which is a good basis to generate a model. Though the value for R 2 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.

Carbon Isotope Ratios-13 C/ 12 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 2 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 2 ) 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 2 , and NO 2 , 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 2 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 2 NO 2 in New York [39,40]. Only China has typically greater NO 2 concentration; in July 2017 the AURA satellite measured 178 billion molecules of NO 2 /mm 2 [39]. Carbon monoxide has its greatest relationship with 13 C/ 12 C in oak in October (R 2 = 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 2 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 2 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 2 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 2 , 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 2 concentration and the carbon isotope ratios of oak was observed (Figures 6 and 10). This may be due to the low resolution of the grid NASA has provided for CO 2 concentration (2 0 × 2.5 0 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 0 × 2.5 0 .
Forests 2020, 11, x FOR PEER REVIEW 17 of 23 Increased CO2 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 CO2, 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 CO2 concentration and the carbon isotope ratios of oak was observed (Figures 6 and 10). This may be due to the low resolution of the grid NASA has provided for CO2 concentration (2 0 × 2.5 0 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 0 × 2.5 0 .

Sulfur Isotope Ratios-34 S/ 32 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]. NO2 concentration is related to pollution [48], March and November NO2 concentrations appear to have some relationship with sulfur isotope ratios with R 2 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 2 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.

Sulfur Isotope Ratios-34 S/ 32 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 2 concentration is related to pollution [48], March and November NO 2 concentrations appear to have some relationship with sulfur isotope ratios with R 2 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 2 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. Figure 11. R 2 of the relationships between sulfur isotope ratios and various monthly mean values for parameters.
The relationship of sulfur isotope ratios, sulfate and precipitation may be due to deposition of sulfate (SO4 2− ) 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 SO4 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 2 between tropospheric NO2 concentration and sulfur ratios in the oak samples. Additionally, when values for NO2 are included in multiple regressions with values for SO4 or RSWR, NO2 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. The relationship of sulfur isotope ratios, sulfate and precipitation may be due to deposition of sulfate (SO 4 2− ) 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 4 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 2 between tropospheric NO 2 concentration and sulfur ratios in the oak samples. Additionally, when values for NO 2 are included in multiple regressions with values for SO 4 or RSWR, NO 2 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.

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 13 C/ 12 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 13 C/ 12 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. Tables 3 and 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 13 C/ 12 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 13 C/ 12 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 13 C/ 12 C isotope ratios across multiple growth rings may be a more practical method despite its apparent limitations.

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.