A δ2H Isoscape of blackberry as an example application for determining the geographic origins of plant materials in New Zealand

Summary

In this investigation, two previously reported precipitation δ²H isoscapes for New Zealand were used to develop a δ²H isoscape for blackberry (Rubus sp.) leaf. These isoscapes were calibrated using the measured δ²H values of 120 authentic blackberry leaf samples collected from across the country. A regression model based on environmental variables available for New Zealand was also determined to predict δ²H values measured from blackberry leaves without initially modelling the precipitation δ²H values. The three models were compared for their accuracy and precision when assigning 10 samples of blackberry leaves for their geographic location based on their measured δ²H values. One of the models based on a precipitation isoscape was similar in accuracy and precision of assignment to the model determined from the environmental variables and provides an approach for determining valid isoscapes for future plant materials.

Keywords:

Leaves – Islands – Environmental geography – New Zealand – Meteorology – Rain – Natural materials

Introduction

Geographic source determination of biological materials based on water stable isotopes has been widely applied in the fields of ecology, anthropology and forensic science [1]. These isotope-based geographic assignment approaches are based upon querying a spatial model of the systematic changes of isotope ratios across a region, continent or the world. Such spatial isotopic models, or isoscapes as originally termed by West et al. [2], represent the changing relative abundances of stable isotopes (e.g. δ²H and δ¹⁸O values) due to changing environmental factors, such as differences in climatological, hydrological or ecological conditions.

The basic requirement for creating an isoscape is an isotopic dataset which represents the spatial area of interest. The Global Network for Isotopes in Precipitation (GNIP) was started in 1961 by the International Atomic Energy Agency in cooperation with the World Meteorological Organization (IAEA-WMO) [3,4]. GNIP provided spatial and temporal information regarding the distribution of ²H, ¹⁸O and ³H in precipitation from more than 1000 stations worldwide, among which 300 stations are engaged for stable isotopes [5]. This data set has allowed for the development of global isoscapes for the H and O isotopic composition of precipitation based on various empirical and geostatistical models [6–10].

While GNIP provides information to estimate global isotopic trends in precipitation, only two GNIP stations have been active in New Zealand [4,6]. These stations were active at Kaitaia (northern end of the North Island) and Invercargill (southern end of the South Island) from 1962–1994 and 1961–2015 respectively [4]. The information from these two stations was insufficient for the purpose of determining localised changes in the isotopic composition of precipitation within the country due to New Zealand’s variable climatic conditions [11]. The Cross-Departmental Research Project (CDRP) [12], funded by the New Zealand Ministry of Agriculture and Forestry (MAF), the New Zealand Department of Conservation (DoC) and the New Zealand Combined Law Agency Group (CLAG) was undertaken to measure a higher spatial and temporal resolution of δ²H values and δ¹⁸O values in precipitation across the country [12]. The CDRP project involved the collection of monthly precipitation from 58 sites in New Zealand and measurement of δ²H and δ¹⁸O values during the period August 2007 ‒ December 2009 inclusive [11–13]. The data obtained through the CDRP project have since been utilised to determine local precipitation isoscapes across New Zealand. Rogers et al. [14] reported an initial mean annual precipitation δ²H isoscape which was calculated using a multiple linear regression (MLR) approach based on elevation, mean annual temperature, and mean annual precipitation. The regression model reported by Rogers et al. [14] was shown to be significant with r² = 0.86 and model parameter estimates were significant at p < α = 0.01. The determined precipitation isoscape was then linked to δ²H values measured from historic feathers from known locations by simple linear regression (r² = 0.48) and utilised to investigate the potential geographic origins of feathers from historic Maori cloaks [14].

It has been shown that temporal variation in environmental parameters like rainfall, temperature and relative humidity, as well as geographic variation in temperature related to latitude and elevation, are drivers of the isotopic composition of precipitation [8]. The addition of measured environmental parameters as ancillary variables into isoscape models has shown to provide improved estimates of isotopic composition in precipitation [8] beyond using purely geostatistical approaches. Baisden et al. [11] reported a revised isoscape for δ²H values in precipitation across New Zealand based on the monthly precipitation δ²H values from the CDRP network, two geographic variables (latitude, elevation) and five daily climate variables (mean sea level pressure, minimum temperature, soil temperature at 10 cm, radiation and wind speed). The monthly weighted means of the daily climate variables were regressed against all available monthly δ²H values from the CDRP network. For monthly predictions, the model reported by Baisden et al. [11] yielded r² = 0.43, for long term annual average predictions the model yielded r² = 0.79.

The precipitation δ²H model reported by Baisden et al. [11] was used by Ehtesham et al. [13] to calibrate via linear regression to a model for determining the δ²H values of milk powders produced in New Zealand. The precipitation δ²H values predicted from the Baisden et al. [11] model were shown to be strongly related to the bulk δ²H value of whole milk powder (r² = 0.86) as well as the compound-specific δ²H values of C14:0, C16:0 and C18:1 fatty acids extracted from milk powder. The precipitation model reported by Baisden et al. [11] allowed for making statistically valid predictions of isotopic composition in precipitation across specific time periods based on actual climate data [11].

Previous work using isoscapes for the geographic assignment of biological materials has shown the need for accurate precipitation isoscapes to underpin the semi-parametric Bayesian framework, which is widely used to determine the likely geographic origins [15]. Calibration of the precipitation isoscape is typically undertaken by determining functions which relate the precipitation isotope values to those of the biological material of interest [15–17]. These calibration functions are usually determined by the evaluation of empirical data collected on known origin samples of the material [15]. The δ²H values of cellulose derived from terrestrial plants have been shown to have a linear relationship with δ²H values in precipitation and can be very simply applied to calibrate a precipitation δ²H isoscape [1,18–25]. However, the relationship between the δ²H values of cellulose in plants and precipitation can be dependent on physical factors like evapotranspiration, temperature and humidity, and as such, some semi-mechanistic models have also been applied [24,25].

The utility of hydrogen isotope ratios as geographic origin indicators has been shown in numerous applications, including archaeology [26], food authentication [27] and forensics [24,28–31]. Booth et al. [29] utilised δ²H as one of several stable isotope indicators for the geographic origins of marijuana grown within and from outside of Alaska. While models based on precipitation δ²H values were determined by Hurley et al. [30,32] and utilised to indicate the origins of marijuana leaf in the conterminous United States. New Zealand specific examples include the previously mentioned application to bird feathers from historical Maori cloaks [14] and to milk and milk powders [13]; the δ²H values in beer and cider samples from New Zealand have been found to reflect the source water of precipitation and irrigation [33,34]. In studies by Holder et al. [35,36], the δ²H values of insect tissues were applied successfully to find information regarding their origin in New Zealand. Development of New Zealand specific isoscapes for plant materials could be utilised for forensic or biosecurity applications, including providing intelligence about growing regions in New Zealand or distinguishing imported plant materials. However, the limitation to this is the ability to obtain authentic samples of plant material with known geographical coordinates with which to calibrate a precipitation isoscape. Ideally, to calibrate the precipitation isoscape, species-specific isotope data from many plant samples covering the geographical extent of interest and inhabiting isotopically different locations would be obtained [16,37].

In this study, the δ²H values measured for blackberry (Rubus sp.) leaf samples from 130 different sites in New Zealand were collected. These samples were used in the development of a blackberry leaf isoscape as a function of δ²H values in precipitation. Blackberry is a mainly uncultivated species introduced to New Zealand that is prevalent across the whole country. It is deciduous, shedding its leaves every winter and is considered a moderately serious pest species in New Zealand as it is threatening native forests and wetlands [38]. This means blackberry is easily found and accessible for sampling across the country, allowing for a survey of plant material from a diverse number of environments.

The overall aim of this investigation was to determine an approach for the development of a plant material isoscape specific to New Zealand, which can be utilised for geographic assignment purposes. The precipitation isoscape models reported by Rogers et al. [14] and Baisden et al. [11] are used to provide the precipitation δ²H values for development of an isoscape which can be used to estimate blackberry leaf δ²H values. The resulting δ²H isoscapes for blackberry leaf are validated and compared to determine the more applicable base precipitation δ²H model. A blackberry leaf isoscape was also determined without prior modelling of the precipitation δ²H values by regression of the measured blackberry leaf δ²H values directly against environmental and geographic variables for further comparison. Prediction of the regions of likely origin for 10 blackberry leaf samples not previously utilised in the modelling is undertaken to provide an example of potential application. The approach taken could then be extended to other plant materials of interest, like marijuana for forensic purposes, or apples and kiwifruit as plant materials of economic and biosecurity importance to New Zealand.

Methods

Sample collection and preparation

A spatially representative set of blackberry samples (n = 130) from the North and South islands of New Zealand were obtained for the study over the period from October 2012 to January 2013. Approximately 30 g of each blackberry leaf sample was picked by hand and collected in a ziplock bag. The blackberry leaf samples were collected from recent, green and uniformly large leaves growing on the top of mature plants. Three individual plants were sampled at each site, and their geographic coordinates were recorded in decimal degrees (WGS84) using a mobile GPS unit (Garmin GPSmap 60CSx). The blackberry collection sites (Fig 1) were situated in areas under the authority of New Zealand local government or the New Zealand DoC. No specific permissions were required for access to these sites or for sample collection as these sites are open to public access and blackberry is not a threatened or protected species in New Zealand.

The samples were transported chilled to the laboratory where possible and stored refrigerated to prevent decomposition until they were able to be dried. The blackberry leaf samples were oven-dried at 30°C for 24 hours [39], which reduced the samples to ~68% of their original mass. The blackberry leaf samples were then stored at -20°C in ziplock bags within black polythene bags until further preparation was required. Approximately 1 g of blackberry leaf was sub-sampled into a MM 400 ball mill (Retsch, Haan Germany) container. The containers were milled overnight resulting in a fine flour. The samples were sieved (250 μm) to obtain a homogenous powder and stored in 10 mL sample vials in a desiccator ready for analysis.

δ²H measurement

The δ²H value of the bulk leaf material was determined using two-stage equilibration [40,41]. The ambient temperature method with vacuum drying was used instead of the steam equilibration method in the analysis of the non-exchangeable δ²H value of the bulk leaf samples. This is because it has been reported that equilibration at ambient temperature produced more consistent and realistic non-exchangeable δ²H values in samples compared to the steam equilibration method [41].

Two sets of approximately 0.6 mg samples in triplicate were weighed into silver capsules (OEA Laboratories, Cornwall UK). A set of triplicates was loaded into each of two separate equilibration chambers along with triplicates of three quality control materials. The silver capsules were closed yet not crimped to allow the exchange of water vapour and effective drying prior to measurement. The samples were equilibrated in the separate equilibration chambers for six days with ~2 mL of different isotopic reference waters (High δ²H value water = +60 ‰ _VSMOW and Low δ²H value water = -262 ‰ _VSMOW) at ~20°C after evacuating the desiccators for 2 min with a vacuum pump to remove air and ensure high concentration of water vapour in the chamber atmosphere.

After equilibration (>6 days), the samples were transferred to autosampler carousels that were placed inside aluminium chambers to isolate from the atmosphere. The samples were dried at 60°C for four days with online evacuation through a freeze-dryer [41]. All the samples were analyzed using a Delta V IRMS (Thermofisher Scientific) with attached TC/EA and Costech zero-blank autosampler. Measured values were reported vs VSMOW. The samples were standardized to reference materials USGS53 (δ²H_vsmow = +40.2 ‰), USGS47(δ²H_vsmow = -150.2± 0.5 ‰) and IAEA-CH-7 (δ²H_vsmow = -100.3 ‰). The quality control was conducted by applying three reference materials, Kudu Horn (δ²H_vsmow = -34.9 ± 0.6 ‰), Caribou Hoof (δ²Hvsmow = -156.4 ± 1.8 ‰) and RH-B (chicken feather) (δ²H_vsmow = 50 ± 2.8 ‰). Note: the original values for these keratin QC standards have been corrected according to the equation presented by Soto et al. [42]. The values presented here are our measurements. A set of the reference materials were measured at the beginning and at the end of every sample run. Duplicate IAEA-CH-7 aliquots were run after every 12 sample aliquots. Note that the QC materials used are all keratinaceous and not plant material. These were the only natural materials available at the time that had internationally accepted values for their non-exchangeable δ²H. As these materials have a very different matrix to the plant samples the results cannot be used for data correction, they are used to monitor the equilibration experiments and to show that consistent results are being obtained in each experiment. The lack of internationally accepted standards of the appropriate matrix and isotope value is a major problem for stable isotope analysis and one that needs to be addressed by the international community. The measured δ²H values are for the 130 blackberry leaf samples are reported in S1 Table.

Statistical modelling

The resultant data set was split into a training set (n = 120) to investigate the variation of the δ²H values as a result of different climatic and environmental conditions in New Zealand. The remaining samples (n = 10) were set aside to be used as an external test set for assignment. The 10 test samples were selected using a stratified random sampling approach. The blackberry samples were grouped by region within New Zealand (defined by New Zealand Regional Council boundaries available from https://datafinder.stats.govt.nz/layer/92204-regional-council-2018-generalised/) and randomly selected from the different groups whilst ensuring no regional group was sampled more than once, to a maximum of 10 selected samples. This ensured the external test set included a spread of different locations and environments across New Zealand to be tested. The external test set was also chosen to represent a scenario similar to that of investigations in food authentication or forensics, where only a small number of independent suspect samples may be available for testing and are in need of being distinguished from particularly identified areas.

Two precipitation δ²H isoscape models have been published for New Zealand. The first model was reported by Rogers et al. [14], the second model was reported by Baisden et al. [11]. The reported underlying datasets for these models (CDRP [12], WorldClim [43] and VCSN [44]) were obtained, and the models calculated in line with each reported method to accurately reproduce each of the precipitation isoscapes as raster grids.

The CDRP data has been previously reported by Frew et al [12] and is available in S2 Table.

The precipitation δ²H isoscape model reported by Rogers et al. [14] utilised data for elevation, mean annual temperature and mean annual precipitation obtained from the WorldClim database as predictor variables. The WorldClim database currently provides open access to 19 bioclimatic variables as an average for 1970–2000 at 10, 5, 2.5 or 0.5 arc minutes spatial resolution for global land areas [43,45]. Rogers et al. [14] utilised an older version of the WorldClim database (version not reported) which was an average of bioclimatic variables over the period 1950–2000 with a resolution of 2.5 arc minutes [46]. The mean annual precipitation data were transformed using the base 10 logarithm before modelling.

The precipitation δ²H isoscape model reported by Baisden et al. [11] was based on multiple linear regression of 2 geographic variables (latitude and elevation), and 5 daily precipitation weighted climate variables (mean sea level pressure, soil temperature at 10 cm, radiation, minimum temperature and wind speed). The climate variables were determined from modelled daily climate data obtained from the Virtual Climate Station Network (VCSN). The VCSN is a climate database, developed and maintained by the New Zealand National Institute of Water and Atmospheric Research (NIWA) [44], that uses interpolation techniques to create a daily climate record at every point on a ~5 × 5 km grid covering all of New Zealand [47,48]. Daily climate variables were weighted by daily precipitation amounts to generate monthly weighted climate means averaged to cover the years 2007–2013, this includes the years of CDRP collection (2007–2009) through to the years of blackberry collection (2012–2013).

Each of these two isoscape models was utilised to determine the δ²H values for precipitation at the blackberry leaf collection sites, utilising an approach similar to that reported by Ehtesham et al. [13]. The average δ²H values in precipitation for grid cells within a 5 km radius of each blackberry sampling site were regressed against the measured δ²H values in the blackberry leaves. Each of the resulting simple linear regression equations, Model 1, which used values from the precipitation δ²H isoscape model of Rogers et al. [14] and Model 2, which used values from the precipitation δ²H isoscape model of Baisden et al. [11], were used to calibrate the precipitation isoscapes for estimation of δ²H values in blackberry leaves.

Calibration of precipitation isocapes to a biological material of interest is mainly undertaken due to a lack of the biological material δ²H values to model directly [1]. Given a reasonable number of biological material samples for this investigation (n = 120) an attempt at modelling the δ²H values of the blackberry leaves against the geographical variables of latitude, longitude (decimal degrees) and altitude (metres above sea level), as well as daily precipitation-weighted environmental variables from the VCSN database by multiple linear regression was also undertaken for comparison purposes (Model 3). Backwards stepwise selection of variables was undertaken following the method outlined in Baisden et al. [11] and the model resulting in the lowest Akaike Information Criterion (AIC) value was selected for use, taking into account the significance of variables entered into the model.

Posterior probability surfaces for geographical assignment were determined from each of the three different models using a semi-parametric Bayesian framework as has been reported in a number of studies and used in several platforms that have been developed for the determination of isoscapes, like IsoMAP and IsoriX [15,49]. In short, Bayes theorem is applied to evaluate the probability, A_i that a sample is from site i given that sample’s measured isotopic value, δ_s, as shown in Eq 1.

Where P(A_i) is the prior probability associated with location i and, the summation is across all locations [15]. For this study, P(A_i) was assumed uniform across all locations. P(δ_s|A_i) is evaluated as shown in Eq 2.

Where δ^i is the isoscape predicted value for the blackberry δ²H at location i and σ_i² is the total variance on the predicted sample value for site i. A similar approach to that outlined by Courtiol et al [49] was followed to determine σ_i², where this total variance depends upon (i) the prediction variance associated with the precipitation model, which varies spatially, (ii) the residual variance of the calibration, which is assumed constant, (iii) the prediction variance of the calibration model (which decreases with the size of the calibration dataset) given the predicted isoscape value at the candidate location and (iv) a covariance term between the prediction error of the calibration model given the predicted isotope value at the candidate location and the prediction error of the precipitation isoscape at this location.

Initial testing of the geographic assignment was undertaken using a repeated k-folds cross-validation approach. For each of 10 repeats, the 120 blackberry samples of the training set were split into 10 folds using random sampling without replacement. Each of the k folds in the repeat was used as a testing set while the remaining k-1 folds were used to redetermine the linear regression function. This allowed for 10 times 10 different folds to be tested (1200 individual tests) whilst ensuring every sample was tested an equal amount of times. A k = 10 fold approach was utilised as this has been empirically shown to limit bias of the test estimate introduced by removal of data from the model building process while incorporating a modest amount of variance [50]. The exact same determined folds were utilised to cross-validate and determine posterior probability surfaces from each of the three models in the investigation. The models were then used to determine the geographical assignments of the external test set of ten samples that had not been utilised in the model building step.

The determined posterior probability surfaces from the cross-validation and external test set were normalised so that all sample probabilities summed to 1 across the spatial domain. The surfaces were then rescaled between 0 and 1 for clarity and ease of comparison, similar to Wunder et al. [37] and Vander Zanden et al. [51]. The resulting surface consists of a range of probabilities where a value close to 0 is a low probability for geographic assignment of the test sample, and a value close to 1 is a high probability for geographic assignment. Areas with rescaled probability >0.667 were considered the likely origin of the test sample, based on the commonly used 2:1 odds ratio threshold reported in previous investigations [17,52].

The efficacy of geographic assignment for the different blackberry leaf δ²H isoscapes were compared using accuracy, precision and similarity metrics as described in Vander Zanden et al. [51]. Individual-level accuracy was evaluated as the relative change in probability between assignments at the known origin site and reflects whether the relative likelihood of origin increases or decreases at the true location when comparing one isoscape to another. Individual-level precision was evaluated as the change in the area of the posterior probability surface at a relative probability density equal to that at the known origin. This metric demonstrates whether the precision increases or decreases for an individual test sample when comparing one isoscape to another. Two-sided paired t-tests were used by Vander Zanden et al. [51] to determine if the mean individual level accuracy and precision values differ from zero. The similarity index reported by Vander Zanden et al. [51], calculates the number of shared cells between two assignment rasters to the total number of possible cells. A function of relative probability against similarity is determined for each individual which is then integrated to find the area under the curve as a value between 0 (no parts match) and 1 (all parts match exactly). This similarity index can be informally considered as the proportion of overlap of the two assignments over the entire range of relative probabilities [51].

All statistical modelling and assessments were undertaken using the R project for statistical computing [53].The R scripts used to conduct assignments and calculate accuracy, precision and similarity have been previously reported in Vander Zanden et al. [51]. A list of the R packages used in this investigation is provided in S3 Table.

Results and discussion

Model 1

The resulting parameters for Model 1, the simple linear regression of measured blackberry δ²H values against precipitation δ²H values using the model reported by Rogers et al. [14] are shown in Table 1.

A plot of the regression result, measured against predicted δ²H values of blackberry leaf exhibits divergent non-linear behaviour at higher predicted values (Fig 2).

This behaviour results in a root mean squared error of 8.6 ‰ and r² of 0.49. Upon further investigation, the non-linear behaviour observed is a result of an approximately 10 ‒ 20 ‰ difference between the measured δ²H values of blackberry from the West Coast of the South Island (as shown by the open circles in Fig 2) and the predicted values for the same region derived from Model 1. The isoscape of model-derived predictions for blackberry leaf δ²H values based on Model 1 are shown in Fig 3.

It can be observed in Fig 3 that along the West Coast of the South Island and across the northern half of the North Island the predicted values do not agree well with the measured values, this is due to the difference observed in the underlying data for the West Coast of the South Island which exerts undue leverage on the regression and skews the prediction surface. The cause of the divergent behaviour observed in Fig 2 is surmised to be due to the temporal extent of the WorldClim data that was averaged into the precipitation δ²H model reported by Rogers et al. [14] which covers records from 1950–2000. Historically, a high amount of average annual precipitation has been recorded for the West Coast of the South Island (4000–5000 mm) over this period compared to the rest of New Zealand (1000–2000 mm). The result of using the linear calibration with the longer-term averages of the precipitation δ²H model reported by Rogers et al. [14] gives higher predicted δ²H values for blackberry leaves along the West Coast of the South Island compared to the measured values and lower predicted δ²H values for blackberry leaves from the upper North Island. A map of the residuals for the calibration model is shown in S1 Fig, which highlights the underestimation and overestimation of δ²H values in different areas of New Zealand. Assessment using a test for Moran’s I was significant at p < 0.01 confirming the presence of the observed spatial autocorrelation in the residuals. For the purposes of geographic assignment, the prediction variances were determined by cross-validation of the Roger’s precipitation δ²H model and Model 1 respectively. The residual variance of Model 1 was determined as 73.9 (s.d. of 8.6 ‰) during the model calculation, and as the precipitation model was able to be reproduced, the covariance terms were also able to be determined. The resulting map of variance for each grid cell ranges from 74 to 85 (s.d. of 8.6 ‰ to 9.2 ‰).

The repeated k-folds cross-validation for Model 1 resulted in, at the selected threshold level, 63% of true locations within the determined region of origin. This came with a mean precision of 74 351 km².

The likely area of origin for each of the test set of 10 blackberry samples was determined using the assignment framework previously described. The figures showing the areas of likely origin are shown in S1 File. Of the 10 test samples, two of these samples (1, and 4) have a low determined probability at their actual location of origin. These test samples are located in the upper half of the North Island of New Zealand, one of the areas affected by the skewed linear calibration. The measured δ²H values for the test samples and those predicted by Model 1 are shown in Table 2 along with the difference between the two δ²H values (measured–predicted) as well as the scaled probability from the geographic assignment.

This result reinforces the importance of having an appropriate precipitation isoscape model, which is related in process to the sample material to be investigated. In this case, due to the higher δ²H values estimated for the precipitation on the West Coast of the South Island and the leverage exerted on the linear regression, the estimates of Model 1 are skewed and not appropriate for the material being investigated.

Model 2

Given that blackberry is deciduous, and sheds leaves in winter, leaf material will, at most, integrate δ²H values for a finite time period prior to sampling (depending on the time of season sampled). This suggests that a precipitation δ²H model based on long term averages may not be appropriate, and an approach which can provide δ²H values across a shorter time period will be more accurate [17]. This has been hypothesised and investigated previously by Vander Zanden et al. [51]. In this previous investigation, the biological materials of interest were monarch butterflies and Eurasian reed warblers. The resulting comparisons showed little to no improvement using a targeted short term precipitation isoscape when compared to a long term isoscape; however there was uncertainty around the effect of water memory assumptions used to estimate the shortened timeframe for monarch butterflies and acknowledgment that short term isoscapes may result in greater assignment efficacy in regions where greater annual variation of the isotopic composition of precipitation may occur [51]. New Zealand is subject to varying exposure from sub-antarctic and subtropical air masses resulting in a climate which is highly variable, and the subject of a more temporally refined isoscape precipitation model put forward by Baisden et al. [11]. The resulting parameters for Model 2, the simple linear regression of measured blackberry leaf δ²H values against precipitation δ²H values using the model reported by Baisden et al. [11] are shown in Table 3.

The non-linear behaviour observed previously in Fig 2 is not present in the plot of the regression result for Model 2 (Fig 4).

This is possibly due to the precipitation model representing a more appropriate time frame for the incorporation of the precipitation δ²H into the sample material. This would suggest that a model based on a shorter timeframe is likely to be more representative of the integration of δ²H values into this type of material as it will not average in the large variances in climatic conditions experienced by New Zealand over longer periods.

The isoscape of model-derived predictions for blackberry leaf δ²H values based on Model 2 is shown in Fig 5.

It can be observed in Fig 5 that the predicted values empirically show better agreement with the measured values than was observed in Fig 3, especially along the West Coast of the South Island and across the northern half of the North Island. This observation is further supported by the smaller RMSE of 4.0 ‰ when the predicted results are compared to the measured results for Model 2 compared to the RMSE of 8.57 ‰ for Model 1. A map of the residuals for Model 2 is presented in S2 Fig, which shows a more random spread across the country than the residuals of Model 1. Assessment using a test for Moran’s I was not significant with p = 0.29, providing little evidence to reject the null hypothesis for the presence of spatial randomness in the residuals. For the geographic assignments, the prediction variances were determined by cross-validation of the Baisden precipitation model and calibration model, respectively. The residual variance of Model 2 was determined as 14.38 (s.d of 3.8 ‰) during the model calculation, and as the precipitation δ²H model was able to be reproduced, the covariance terms were also able to be determined. The resulting map of variance for each grid cell ranges from 14 to 24 (s.d. of 3.7 to 4.9 ‰).

The repeated k-folds cross-validation for Model 2 resulted in, at the selected threshold level, 68% of true locations within the determined region of origin. This came with a mean precision of 31 002 km².

The likely area of origin for each of the test set of 10 blackberry samples was also determined for this model. The figures showing the areas of likely origin are shown in S2 File. 9 of the 10 test samples have a high determined probability at their actual location of origin with the true location of test sample 8, having a probability of 0.56. The measured δ²H values for the test samples and those predicted by the model based on the precipitation isoscape from Baisden et al. [11] are shown in Table 4 along with the difference between the two δ²H values (measured–predicted) as well as the scaled probability from the geographic assignment.

Model 3

Another approach, where sufficient data are available, is to directly model the changing patterns in the target biological material. There are a number of modelling approaches reported in the literature that could be applied to do this, including empirical comparisons using linear regression [8], geostatistical modelling of the semivariance [28] or mechanistic modelling [24]. For comparison purposes, in this investigation, the measured δ²H values for the blackberry leaves were empirically modelled by multiple linear regression (Model 3). The parameters of the resulting MLR model for measured blackberry leaf δ²H values against daily precipitation weighted environmental variables are shown in Table 5.

A plot of the regression result from Model 3 against the measured blackberry leaf δ²H values is shown in Fig 6.

In this case, when comparing the δ²H values predicted by the model with the measured values, the resultant RMSE (3.7 ‰) is similar to that obtained for Model 2 (4.0 ‰). The isoscape of model-derived predictions for blackberry leaf δ²H values based on Model 3 is shown in Fig 7.

A map of the residuals for the model is shown in S3 Fig, which exhibits a random spread across the country. Assessment using a test for Moran’s I was not significant with p = 0.29, similar to that of Model 2, providing little evidence to reject the null hypothesis for the presence of spatial randomness in the residuals. The obvious advantage of empirically modelling the δ²H values of the sample material directly is the reduced sources of variance due to removing the need for modelling precipitation δ²H and then undertaking further regression. The total variance term used for the geographic assignment of the test set was a summation of the prediction variance of Model 3 determined by cross-validation, which varies spatially, and the residual variance of Model 3. This was determined as 13.5 (s.d. of 3.67 ‰). Given that the predictions were not based on the prediction of precipitation δ²H values, there was no need to include the prediction variance and covariance terms associated with this extra modelling step.

The repeated k-folds cross-validation for Model 1 resulted in, at the selected threshold level, 68% of true locations within the determined region of origin. This came with a mean precision of 29 233 km².

The figures showing the likely origin of the test samples are shown in S3 File. 8 out of 10 test samples have a high determined probability at their actual location of origin. The measured δ²H values for the test samples and those predicted by Model 3 are shown in Table 6 along with the difference between the two values (measured–predicted) as well as the scaled probability from geographic assignment.

Model comparisons

Model 3 was compared to the two models (Models 1 and 2), which were calibrated against the different precipitation isoscapes. Validation of the different blackberry leaf δ²HHHH isoscape models was undertaken using repeated k-folds cross-validation as well as the external test set (n = 10), which was segregated prior to the development of the models with the training set data. A comparison of the models was undertaken in terms of their individual accuracy and precision utilising the metrics outlined by Vander Zanden et al. [51].

When comparing the cross-validation results of Models 1 and 2, the individual level assignment accuracy was observed to be slightly lower on average for Model 1. A two-sided paired t-test showed the observed difference in accuracy was not significant with p = 0.03 > α = 0.01 (mean = -0.02, t = -2.13, df = 1 199). A similar result was obtained when comparing Models 1 and 3, the accuracy of Model 1 was observed to be slightly lower however this was shown to not be significant with p = 0.35 > α = 0.01 (mean = -0.01, t = -0.93, df = 1 199). Very little difference in accuracy was determined comparing Models 2 and 3 with p = 0.05 > α = 0.01 (mean = -0.00002, t = 1.9243, df = 1 199). This showed that overall there was no real difference in accuracy across the three models.

The individual-level assignment precision was also undertaken through assessment of the change in assignment area from one model to another in thousands of km². Comparison of Models 1 and 2 showed an average reduction in assignment area of 43 349 km² with Model 2 being the more precise (p < 0.01, t = 22.78, df = 1 199). Comparison of Models 1 and 3 gave a similar result with Model 3 being more precise with an average reduction in assignment area of 45 119 km² (p <0.01, t = 22.41, df = 1 199) over Model 1. Comparison of Models 2 and 3 showed a smaller change in precision with Model 3 having an average reduction in assignment area of only 1 769 km² (p < 0.01, t = 3.72, df = 1 199). While this is still a significant change in precision based on the two-sided paired t-test, the smaller average change in assignment area shows that Models 2 and 3 have a more similar precision to each other than Model 1 has to either of the other two models.

The similarity indices for the cross-validation assignments of Model 1 compared to Model 2 ranged from 0.2 proportion of overlap to 0.6, with a mean proportion of overlap of 0.46. For Model 1 compared to Model 3, the similarity indices had a similar range from 0.2 proportion of overlap to 0.55 with a mean proportion of overlap of 0.41. Given that there was no significant difference in accuracy between Model 3 and Model 2 and a small though significant change in precision, it is expected that Model 2 and Model 3 will have similar proportions of overlap in assignment area. This is reflected in the similarity indices of the comparison of Model 2 to Model 3 which ranged from 0.6 proportion of overlap to 0.85, with a mean proportion of overlap of 0.76. Histograms of the change in relative probability (accuracy), change in area (precision) and similarity indices for each of the model comparisons are provided in S4 File.

When comparing the geographic assignment outcomes of the external test set, it was observed that for test sample 8, Model 1, makes predictions which are closer to the measured value than the predictions of Model 2 and Model 3. Test Sample 8 surprisingly originates from the West Coast area of the South Island of New Zealand for which Model 1 did not behave linearly. Further investigation of this test sample shows it exhibits a higher δ²H value (-75.0 ‰) than the surrounding training samples measured from the same region (~-85 –-95 ‰). While Model 1 gives a better estimate of the mean δ²H value at the location of test sample 8 the efficacy of the geographic assignment is more relevant to food authentication or forensic application as not only the accuracy of the estimation and assignment is important but the precision also. If the intended purpose of the model is to provide a likely area of origin for further investigation, then applying resources across a smaller area is more efficient.

The results of the comparison using the external test set are shown in Table 7; the precision difference is shown as the area difference in thousands of km².

While individual differences in accuracy can be seen across the external test sample assignments, the different models show no significant change in mean accuracy (p-values >0.1) when they are compared, as observed from cross-validation. However, the mean precision of the assignments made by the isoscape based on Model 1 is significantly worse than the precision of the assignments made by the isoscapes based on Model 2 (p = 0.01) and Model 3 (p = 0.02). The difference in the precision of assignments made between the isoscapes based on Model 2 and Model 3 was not significant (p = 0.61) with Model 3 having an average 1 631 km² smaller assignment area than Model 2 across all test samples. This is in agreement with the assessments of the accuracy and precision determined by repeated k-folds cross-validation.

Table 8 shows the calculated similarity indices for the 10 test samples between each of the considered isoscapes.

The similarity indices show the assignment areas from the isoscape based on Model 1 for each test sample are not very similar to those based on Model 2 or Model 3; with mean proportions of overlap of 0.44 and 0.42 respectively. Model 2 and Model 3 show a higher degree of similarity with a mean proportion of overlap of 0.84. These outcomes are in line with those from the repeated k-folds cross validation.

Conclusions

Two previously reported precipitation isoscapes were calibrated to the δ²H measurements of 120 samples of blackberry leaves from across New Zealand (Models 1 and 2). The measured δ²H values of the 120 blackberry leaf samples were also modelled against daily weighted environmental variables (Model 3). The three resulting isoscapes were compared to determine the more appropriate for geographic assignment. Of the three isoscapes, Model 1 which was based on the precipitation isoscape reported by Rogers et al. [14] was the least precise in terms of geographic assignment and exhibited non-linear behaviour in the underlying calibration from precipitation δ²H values to blackberry leaf δ²H values. This resulted in a skewed prediction surface, higher residual error and higher prediction error which gave rise to the lack of relative precision. The non-linear behaviour was surmised to be due to the incorporation of long-term climate information, which was inappropriate for the shorter integration time of δ²H into blackberry leaves. Model 2, which was based on the precipitation isoscape reported by Baisden et al. [11] performed with high accuracy when making assignments of the external test set, including 9 out of 10 true locations within the assignment area. The model also had significantly higher precision when compared to the result of Model 1 based on the undertaken cross-validation and external test set. The precision of Model 2 was not significantly different from the precision of Model 3. The change in precision between Model 2 and Model 3 was considered significant based on the cross-validation result, however. the difference in average area between the two models was approximately 26 times smaller than the average differences between these models and Model 1. Model 3 involved modelling of climate variables to the δ²H values of the blackberry to determine an isoscape for geographic assignment. However, the approach taken to determine Model 3 requires a large number of measurements of the plant material from representative sites across the area of interest to capture adequate variation. This may be difficult for plant materials of a more clandestine nature. Given the agreement in accuracy, precision and the high degree of similarity between Model 3 and Model 2; the approach of calibrating the precipitation isoscape reported by Baisden et al [11] is as viable as modelling plant material δ²H values against the climate parameters. The approach for determining the precipitation δ²H isoscape reported by Baisden et al. [11] is flexible enough that it can represent a range of integration times for the plant material. It appears to be an appropriate choice of underlying model against which specific plant material δ²H values could be calibrated for New Zealand and could be utilised as the starting point to calibrate multiple different plant materials against.

Supporting information

S1 Table [wgs84]
The measured δH values of the 130 blackberry leaf samples.

S2 Table [wgs84]
The New Zealand CDRP isotopes in precipitation data used for determining precipitation isoscapes for New Zealand.

S3 Table [xlsx]
R packages utilised for this investigation and their accompanying citations.

S1 Fig [tif]
Map of the residuals from the linear calibration of δH values in precipitation from Rogers et al. to δH values in blackberry leaf (Model 1).

S2 Fig [tif]
Map of the residuals from the linear calibration of δH values in precipitation from Baisden et al. to δH values in blackberry leaf (Model 2).

S3 Fig [tif]
Map of the residuals from the multiple linear regression of δH values in blackberry leaf to daily precipitation-weighted environmental variables obtained from the NIWA VCSN database (Model 3).

S1 File [black]
Geographic assignment maps for test samples 1 to 10 using the δH blackberry leaf isoscape determined by calibration of the δH precipitation isoscape reported by Rogers et al. (Model 1).

S2 File [black]
Geographic assignment maps for test samples 1 to 10 using the δH blackberry leaf isoscape determined by calibration of the δH precipitation isoscape reported by Baisden et al. (Model 2).

S3 File [black]
Geographic assignment maps for test samples 1 to 10 using the δH blackberry leaf isoscape determined by direct modelling of daily precipitation-weighted environmental variables obtained from the NIWA VCSN database (Model 3).

S4 File [pdf]
Histograms of the change in relative probability (accuracy), the change in area (precision) and similarity indices for each of the model comparisons based on repeated k-folds cross-validation.