Using plant physiological stable oxygen isotope models to counter food fraud
Independent reference samplesThe authentic, independent strawberry (Fragaria × ananassa) reference samples used for model validation in this study were provided by Agroisolab GmbH (Jülich, Germany). The samples were collected either directly by the company or on their behalf through authorized sample collectors between 2007 and 2017. The primary purpose of such authentic reference samples is the direct comparison between their stable isotope compositions (oxygen, hydrogen, carbon, nitrogen, or sulfur) to those of samples of suspect origin. Accompanying metadata for each reference sample included information about the geographic origin as community name, postal code, or location coordinates, and information about the month and year the strawberry sample was picked. In total, we used δ18O values from 154 reference samples. Most samples were collected in the UK, Germany, Sweden, and Finland (Fig. 1). All reference samples were grown on open strawberry-fields rather than artificial greenhouse conditions. All berry samples were collected from cultivated, non-endangered plant species (“garden strawberry”), and the research conducted complies with all relevant institutional, the corresponding national, and also international guidelines and legislation.After collection in the field, samples were stored in airtight containers and shipped directly to Agroisolab, where they were stored frozen prior to analysis. In order to analyze the oxygen stable isotope composition of the organic strawberry tissue, the lipids were solvent-extracted with dichloromethane for a at least 4 h, using a Soxhlet extractor. The remaining samples were dried and milled to a fine powder. 1.5 mg of the powder was weighed into silver capsules. The silver capsules were equilibrated for at least 12 h in a desiccator with a fixed relative humidity of 11.3%. After a further vacuum drying the samples were measured via high-temperature furnace (Hekatech, Wegberg, Germany) in combination with an Isotope-ratio mass spectrometer (IRMS) Horizon (NU Instruments, Wrexham, UK). The pyrolysis temperature was 1530 °C and the pyrolysis tube consisted of covalent-bound SiC (Agroisolab patented). The reproducibility of the measurement was better than 0.6 ‰.Oxygen isotope model calculationPlant physiological stable isotope models simulate the oxygen isotopic composition of leaf water or organic compounds synthesized therein as δ18O values in per mil (‰), where δ18O = (18O/16O)sample/(18O/16O)VSMOW − 1, and VSMOW is Vienna Standard Mean Ocean Water as defined by the VSMOW-Standard Light Antarctic Precipitation (SLAP) scale. The Craig-Gordon model57, which was developed to mathematically describe the isotopic enrichment of standing water bodies during evaporation and later modified for plants, is the basis for modelling plant water δ18O values23,58. Plant source water is the baseline for the model, which is the precipitation-derived soil water that plants take up through their roots without isotope fractionation51,59,60. The 18O enrichment of water within leaves is described by the following equation (Eq. 1)36,61:$$ Delta^{18} {text{O}}_{{text{e leaf}}} = left( {1 +upvarepsilon ^{ + } } right)left[ {left( {1 + {upvarepsilon }_{{text{k}}} } right)left( {1 – {text{e}}_{{text{a}}} /{text{e}}_{{text{i}}} } right) + {text{e}}_{{text{a}}} /{text{e}}_{{text{i}}} (1 + Delta^{18} {text{O}}_{{{text{Vapor}}}} )} right]{-}1 $$
(1)
where Δ18Oe_leaf is the oxygen isotopic enrichment above source of water at the evaporative site in leaves, ε+ is the equilibrium fractionation between liquid water and water vapor, εk is the kinetic fractionation associated with the diffusion through the stomata and the boundary layer. ea/ei is the ratio of ambient vapor pressure in the atmosphere to intercellular vapor pressure in the leaf. Δ18OVapor is the isotopic composition of the ambient vapor above source water, which in this study is assumed to be in equilibrium with the source water (Δ18OV = − ε+)62,63. This assumption can be used, if the atmosphere is well mixed, and plants’ source water derives from recent precipitation events. For crops, growing in the temperate climate of the mid latitudes this is usually the case, especially over the long time periods (several weeks) over which strawberries grow. If such a model is applied in other climatic zones (e.g. tropics), this assumption should, however, be reevaluated64. The equilibrium fractionation factor (ε+)65,66 and kinetic fractionation factor (εk)67 can be calculated with the following equations (Eqs. 2 and 3):$$upvarepsilon ^{ + } = left[ {exp left( {frac{1.137}{{left( {273 + T} right)^{2} }}*10^{3} – frac{0.4156}{{273 – T}} – 2.0667*10^{ – 3} } right) – 1} right]*1000 $$
(2)
where T is the leaf temperature in degrees Celsius. In our calculations, leaf temperature was set to 90% of the monthly mean air temperature, which describes a realistic leaf-energy balance scenario for well-watered crops68,69, and also yielded the best model performance with respect to the reference data. As leaf to air temperature differences have a strong influence on leaf water δ18O values, this assumption needs to be independently tested in future applications. For example, changing leaf temperature from 20 °C to 22 °C at a constant air temperature of 20 °C and a source water δ18O value of -10 ‰ will affect leaf water δ18O values by + 1.4 ‰.$$upvarepsilon _{{text{k}}} = frac{{28{text{r}}_{{text{s}}} + 19{text{r}}_{{text{b}}} }}{{{text{r}}_{{text{s}}} + {text{r}}_{{text{b}}} }} $$
(3)
where rs is the stomatal resistance and rb is the boundary layer resistance in m2s/mmol, which is the inverse of the stomatal and boundary layer conductance. For our model calculations, we consistently used stomatal conductance values of 0.4 mol/m2s, stomatal resistance values of 1 m2s/mol70.The Craig-Gordon model predicted leaf water values are often enriched in 18O relative to measured bulk leaf water δ18O values26,27. This is because the model describes the δ18O values of water at the site of evaporation while measurements typically give bulk leaf water δ18O values36,71. The two-pool modification to the Craig-Gordon model corrects for this effect by separating bulk leaf water into a pool of evaporatively enriched water at the site of evaporation (δ18Oe_leaf derived from the Craig-Gordon model, Eq. 1) and a pool of unenriched plant source water (δ18Osource water)25. δ18Oe_leaf is calculated as follows:$$updelta ^{18} {text{O}}_{{{text{e}}_{text{leaf}}}} = , (Delta^{18} {text{O}}_{{{text{e}}_{text{leaf}}}} +updelta ^{18} {text{O}}_{{text{source water}}} ) , + , (Delta^{18} {text{O}}_{{{text{e}}_{text{leaf}}}} * ,updelta ^{18} {text{O}}_{{text{source water}}} )/1000) $$
(4)
In the two-pool modified Craig-Gordon model (Eq. 5), the proportion of unenriched source water is described as fxylem36.$$updelta ^{18} {text{O}}_{{text{leaf water}}} = , left( {1 , {-}f_{xylem} } right) , * ,updelta ^{18} {text{O}}_{{{text{e}}_{text{leaf}}}} + , left( {f_{xylem} * ,updelta ^{18} {text{O}}_{{text{source water}}} } right) $$
(5)
Values for fxylem in leaf water generally range from 0.10 to 0.3336,37,38,39,40 but higher values have also been observed72. For strawberry plants, leaf water fxylem values were recently shown to vary between 0.24 and 0.3428.Organic molecules in leaves generally reflect the δ18O values of the bulk leaf water plus additional isotopic effects occurring during the assimilation of carbohydrates and post-photosynthetic processes21,22,34. The fractionations occur when carbonyl-group oxygen exchanges with leaf tissue water during the primary assimilation of carbohydrates (trioses and hexoses)42. This process causes 18O enrichment, described as εwc42, and has been determined to be ~ + 27 ‰21,22,73.During the synthesis of cellulose from primary assimilates, sucrose molecules are broken down to glucose and re-joined, allowing some of the carbonyl group oxygen to further exchange with water in the developing cell. The isotopic fractionation (εwc) during this process is assumed to be the same as in the carbonyl oxygen exchange during primary carbohydrate assimilation (~ + 27 ‰)41,42. During the formation of cellulose, the δ18O values of the primary assimilates are thus partially modified by the water in the developing cell33. Equation (6) describes this process34$$updelta ^{18} {text{O}}_{{{text{cellulose}}}} = {text{ p}}_{{text{x}}} {text{p}}_{{{text{ex}}}} *left( {updelta ^{18} {text{O}}_{{text{source water}}} + ,upvarepsilon _{{{text{wc}}}} } right) , + , left( {1 , – {text{p}}_{{text{x}}} {text{p}}_{{{text{ex}}}} } right) , * , left( {updelta ^{18} {text{O}}_{{text{leaf water}}} + ,upvarepsilon _{{{text{wc}}}} } right) $$
(6)
where δ18Ocellulose is the oxygen isotopic composition of cellulose, pex is the fraction of carbonyl oxygen in cellulose that exchanges with the medium water during synthesis, and px is the proportion of unenriched source water in the bulk water of the cell where cellulose is synthesized33. Bulk water in developing cells where cellulose is synthesized, i.e. in the leaf growth-and-differentiation zone, has been found to primarily reflect the isotope composition of source water43. Therefore, px in Eq. 6 is likely larger than fxylem in Eq. (5). For practical reasons, the parameters px or pex are typically not determined individually, but as the combined parameter pxpex45. For cellulose in leaves of grasses, crops, and trees pxpex has been found to range from 0.25 to 0.5445,46,47,48,49.In this study as in many applied examples where plant δ18O values are used for origin analysis we attempt to simulate the δ18O values of dried bulk tissue. Bulk dried plant tissue (δ18Obulk) contains in addition to carbohydrates compounds such as lignin, lipids, and proteins, which can be 18O-depleted compared to carbohydrates50. Since this needs to be accounted for in the model, we included the parameter c into the model. As pxpex and c cannot be determined separately they are used as a combined model parameter in our approach pxpexc:$$updelta ^{18} {text{O}}_{{{text{bulk}}}} = {text{ p}}_{{text{x}}} {text{p}}_{{{text{ex}}}} {text{c}}*left( {updelta ^{18} {text{O}}_{{text{source water}}} + ,upvarepsilon _{{{text{wc}}}} } right) , + , left( {1 – {text{ p}}_{{text{x}}} {text{p}}_{{{text{ex}}}} {text{c}}} right) , * , left( {updelta ^{18} {text{O}}_{{text{leaf water}}} + ,upvarepsilon _{{{text{wc}}}} } right) $$
(7)
Bulk dried tissue δ18O values of strawberries in Cueni et al. (in review) did not differ statistically from pure cellulose δ18O values in strawberries. Consequently, pxpex and pxpexc are identical for strawberries and ranges from 0.41 to 0.51. This approach allows the calculation of bulk dried tissue δ18O values without the knowledge of cellulose δ18O values, which is the case for the data set used in this study, and contrasts to the approach by Barbour & Farquhar (2000), where bulk dried tissue δ18O values are assessed by an offset (εcp) to the cellulose δ18O values.Model parameter selectionTo find the best values of the key model parameters for the prediction of strawberry bulk dried tissue δ18O values, we used different combinations of the values for the parameters. Specifically, we compared average parameter values from the literature that were derived from leaves and parameter values that were specifically derived for berries (Cueni et al. in review) to test if a leaf-level parameterization of the model is sufficient or if a berry-specific parameterization is necessary for producing satisfying model prediction. These values were either (i) fxylem and pxpex values reported in literature for leaf water and cellulose from various species that were averaged, (ii) values averaged for leaves (fxylem) and berries (pxpexc) of berry producing plants, or (iii) values for leaves (fxylem) and berries (pxpexc) specifically obtained for strawberry plants. For the general leaf-derived parameter values we used mean literature values originally obtained for leaf water and leaf cellulose δ18O for different species and averaged these values (0.22 for fxylem36,37,38,39,40 and 0.40 for pxpex45,46,47,48,49) (Table 1). For berries (average of the values of raspberries and strawberries) the mean leaf-derived fxylem value was 0.26 and the value for pxpexc was determined to be 0.46 (Table 1) (data derived from Cueni et al. in review). For strawberry plants, the leaf-derived fxylem value we used was 0.30, and the value for bulk dried tissue (pxpexc) was determined to be 0.46 (Table 1) (data derived from Cueni et al. in review). Since the pxpexc values of different berry species did not differ, this resulted in a total of six different model input parameter combinations.Table 1 Values of the model parameters (fxylem and pxpex/pxpexc) used for the simulations of strawberry bulk dried tissue δ18O values.Full size tableEnvironmental model input data selectionIn order to apply the strawberry parametrized bulk dried tissue oxygen model on a spatial scale, spatially gridded climate and precipitation isotope data layers were used as model inputs. The accurate simulation of geographically distinct δ18O values, however, requires the use of the most appropriate and best available input variables. We therefore tested the importance of the temporal averaging and lead time of the input data relative to the picking date of the berry. We defined these collectively as the “integration time” of the input data. Climate of the growing season46,53, and precipitation δ18O values of rain-events prior and during the growing season51,54,55 have been shown to shape plant tissue water and organic compound δ18O values. The major objective of our study was thus a careful evaluation of the most appropriate type and integration time of model input variables needed for this kind of model simulation. Moreover, to find the best data source provided, we also used several different spatial climate and precipitation isotope datasets in our evaluations (Table 2).Table 2 (a) Table showing the different climatic (air temperature and vapor pressure) and isotope (precipitation and vapor δ18O) data products, (b) as well as the different integration times and names used in the study used to simulate strawberry bulk dried tissue δ18O values.Full size tableTwo precipitation isotope data products were compared (Table 2): (1) The mean monthly precipitation δ18O grids by Bowen (2020), which are updated versions of the grids produced by Bowen and Revenaugh (2003) and Bowen et al. (2005) (Online Isotopes in Precipitation Calculator, OIPC Version 3.2). They provide global grids of monthly long-term mean precipitation isotope values. The resolution of these global grids is 5’. (2) Precipitation isotope predictions from Piso.AI (Version 1.01)32. This source provides values for individual months and years based on station coordinates32. Both data sets were on the one hand used for the precipitation δ18O input data of the model, and also to extrapolate the vapor δ18O values from sets (see model description above), which we treated as two individual, independent input data sets.For the climatic drivers of the model (air temperature and vapor pressure), we used the gridded data products from the Climatic Research Unit (CRU) (TS Version 4.04)29 and the E-OBS gridded dataset by the European Climate Assessment & Dataset (Version 22.0e)30 (Table 2). The CRU dataset provided global gridded monthly mean air temperature, and mean vapor pressure with a resolution of 0.5°. The E-OBS dataset included European daily mean air temperature and relative humidity gridded data, with a resolution of 0.1 arc-degrees. We calculated monthly mean air temperature and relative humidity grid layers, on the basis of these daily mean air temperature and relative humidity grids, respectively.Fruit tissue formation takes place over a period of several weeks leading up to picking date46,53. This results in a lead time between the date that best represents the mean climate conditions, and source water and vapor stable isotope signal influencing the isotope signal during tissue formation, and the picking date. As integration time of the input data, we therefore investigated lead times of 1, 2, 3, and 4 months, as well as the three months leading up to the picking date (Table 2). Moreover we also used more general European strawberry growing season averages76, independent of either the sampling month (yearly May to July mean) or the sampling year (2007 to 2017 May to July mean) (Table 2). Precipitation isotope data means were calculated as amount-weighted averages using CRU mean monthly precipitation data. This means that the long term mean precipitation δ18O values taken from OIPC were weighted by yearly specific CRU monthly precipitation totals for the case of the three months or growing season averages for individual years, and by average monthly precipitation totals (May, June, and July) from 2007 to 2017 for the long-term growing season calculation. The same assessment was also made using precipitation values from Piso.AI.Validation of model with reference samplesUsing the plant physiological model described above, we calculated the strawberry bulk dried tissue δ18O values for the location and the growing time of each authentic reference sample. For the model input data, we tested variable combinations using each of the eight integration times described in Table 2, along with all combinations of the data sources outlined in Table 2. This resulted in a total of 65,536 combinations of input variables per model parameter combination (fxylem and pxpex/pxpexc, Table 1), yielding model results to be evaluated against the measured reference samples. Our approach can be described with the following equation:$$updelta ^{18} {text{O}},{text{plant }} = fleft( {{text{air}},{text{temp}}left( {text{s,t}} right),{text{ relative}},{text{humidity}}left( {text{s,t}} right), updelta ^{18} {text{O}},{text{precip}}.left( {text{s,t}} right),updelta ^{18} {text{O}},{text{vapor }}left( {text{s,t}} right)} right) $$
(8)
where δ18O plant is the simulated δ18O value of the strawberry, s is the data product for the specified input variable (Table 2), and t is the integration time of the specified variable (Table 2).For the crucial model parameters fxylem and pxpex/pxpexc we on the one hand used the values proposed for leaves by literature (for pxpex), and on the other hand average of the values of raspberries and strawberries and strawberry-specific values determined from Cueni et al. (in review) (for pxpexc). In all calculations an εwc value of + 27 ‰ was used. To calculate mean monthly relative humidity values from the provided CRU vapor pressure data, site specific elevation was extracted from the ETOPO1 digital elevation model77, and used to calculate the approximate atmospheric pressure. These values were then used in combination with air temperature to calculate the saturation vapor pressure after Buck (1981), in order to assess relative humidity (relative humidity = vapor pressure/saturation vapor pressure). The R-script of the model is available on “figshare”, find the URL in the data availability statement.Statistical analysesStatistical analyses were done using the statistical package R version 3.5.379. The relationships of the range δ18O values observed with latitude, and between CRU and E-OBS mean air temperature were compared with a linear regression model, and with an alpha level that was set to α = 0.05. The results of the 65,536 models for each of the six physiological parameter combinations were compared with the measured δ18O bulk dried tissue values of the authentic reference samples (n = 154) by calculation of the root mean squared error (RMSE).Calculation of prediction mapsPrediction maps showing the regions of possible origin of a sample with unknown provenance are the product that is of interest in the food forensic industry. We calculated the prediction maps shown in Fig. 5 for three example δ18O values of strawberries collected in July 2017: (i) + 20 ‰ representing a mean Finish/Swedish sample, (ii) + 24.5 ‰ representing a mean German sample, and (iii) + 27 ‰ representing a mean southern European sample.The prediction maps were calculated in a two-step approach. First, we calculated a map of the expected strawberry bulk dried tissue δ18O values of berries grown in July 2017. For this we used the average berry model input parameters (fxylem and pxpexc, Table 1), and the best fitting model input data and integration time combination, which we assessed beforehand (Fig. 2). We thus used CRU mean air temperature and vapor pressure from June 2017, precipitation δ18O values from OIPC as an average from April, May and June, and vapor δ18O values calculated from OIPC precipitation δ18O values from April. Since using spatial maps as model input data, this calculation resulted in a mapped model result. In a second step, we calculated the prediction maps. For this we first subtracted the δ18O value of the bulk dried tissue of the sample strawberry from the mapped result of the best berry-specific model. This was done for each pixel value of the map. This resulted in a map showing the difference of the sample δ18O value and the predicted map δ18O value for each pixel of the map. The places (pixels) that are predicted to have the same δ18O value as the sample strawberry thus are represented by a value of zero. Based on the prediction error of the best berry-specific model (RMSE = 0.96 ‰), the one sigma (68%) and two sigma (95%) confidence intervals around the areas showing no difference to the δ18O value of the suspected sample could be assessed. This means that the bigger the difference between the simulated δ18O value and the sample value, the lower the probability of provenance of the sample. In other words, a difference between the sample’s δ18O value and the predicted δ18O value of 0 ‰ to ± 0.96 ‰ equals a possible provenance of at least 68% (one sigma), and a difference between ± 0.96 ‰ and ± 1.92 ‰ reflects a possible provenance between 68 and 27% (two sigma). Regions on the map with bigger differences than ± 1.92 ‰ represent regions of possible provenance, lower than 5% (bigger than two sigma). More
