Using plant physiological stable oxygen isotope models to counter food fraud

Independent reference samples

The 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 δ¹⁸O 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 calculation

Plant physiological stable isotope models simulate the oxygen isotopic composition of leaf water or organic compounds synthesized therein as δ¹⁸O values in per mil (‰), where δ¹⁸O = (¹⁸O/¹⁶O)_sample/(¹⁸O/¹⁶O)_VSMOW − 1, and VSMOW is Vienna Standard Mean Ocean Water as defined by the VSMOW-Standard Light Antarctic Precipitation (SLAP) scale. The Craig-Gordon model⁵⁷, 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 δ¹⁸O values^23,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 fractionation^51,59,60. The ¹⁸O enrichment of water within leaves is described by the following equation (Eq. 1)^36,61:

where Δ¹⁸O_{e_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. e_a/e_i is the ratio of ambient vapor pressure in the atmosphere to intercellular vapor pressure in the leaf. Δ¹⁸O_Vapor 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 (Δ¹⁸O_V = − ε⁺)^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 reevaluated⁶⁴. The equilibrium fractionation factor (ε⁺)^65,66 and kinetic fractionation factor (ε_k)⁶⁷ can be calculated with the following equations (Eqs. 2 and 3):

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 crops^68,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 δ¹⁸O 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 δ¹⁸O value of -10 ‰ will affect leaf water δ¹⁸O values by + 1.4 ‰.

where r_s is the stomatal resistance and r_b is the boundary layer resistance in m²s/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/m²s, stomatal resistance values of 1 m²s/mol⁷⁰.

The Craig-Gordon model predicted leaf water values are often enriched in ¹⁸O relative to measured bulk leaf water δ¹⁸O values^26,27. This is because the model describes the δ¹⁸O values of water at the site of evaporation while measurements typically give bulk leaf water δ¹⁸O values^36,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 (δ¹⁸O_{e_leaf} derived from the Craig-Gordon model, Eq. 1) and a pool of unenriched plant source water (δ¹⁸O_{source water})²⁵. δ¹⁸O_{e_leaf} is calculated as follows:

In the two-pool modified Craig-Gordon model (Eq. 5), the proportion of unenriched source water is described as f_xylem³⁶_.

Values for f_xylem in leaf water generally range from 0.10 to 0.33^{36,37,38,39,40} but higher values have also been observed⁷². For strawberry plants, leaf water f_xylem values were recently shown to vary between 0.24 and 0.34²⁸.

Organic molecules in leaves generally reflect the δ¹⁸O values of the bulk leaf water plus additional isotopic effects occurring during the assimilation of carbohydrates and post-photosynthetic processes^21,22,34. The fractionations occur when carbonyl-group oxygen exchanges with leaf tissue water during the primary assimilation of carbohydrates (trioses and hexoses)⁴². This process causes ¹⁸O enrichment, described as ε_wc⁴², 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 δ¹⁸O values of the primary assimilates are thus partially modified by the water in the developing cell³³. Equation (6) describes this process³⁴

where δ¹⁸O_cellulose is the oxygen isotopic composition of cellulose, p_ex is the fraction of carbonyl oxygen in cellulose that exchanges with the medium water during synthesis, and p_x is the proportion of unenriched source water in the bulk water of the cell where cellulose is synthesized³³. 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 water⁴³. Therefore, p_x in Eq. 6 is likely larger than f_xylem in Eq. (5). For practical reasons, the parameters p_x or p_ex are typically not determined individually, but as the combined parameter p_xp_ex⁴⁵. For cellulose in leaves of grasses, crops, and trees p_xp_ex has been found to range from 0.25 to 0.54^{45,46,47,48,49}.

In this study as in many applied examples where plant δ¹⁸O values are used for origin analysis we attempt to simulate the δ¹⁸O values of dried bulk tissue. Bulk dried plant tissue (δ¹⁸O_bulk) contains in addition to carbohydrates compounds such as lignin, lipids, and proteins, which can be ¹⁸O-depleted compared to carbohydrates⁵⁰. Since this needs to be accounted for in the model, we included the parameter c into the model. As p_xp_ex and c cannot be determined separately they are used as a combined model parameter in our approach p_xp_exc:

Bulk dried tissue δ¹⁸O values of strawberries in Cueni et al. (in review) did not differ statistically from pure cellulose δ¹⁸O values in strawberries. Consequently, p_xp_ex and p_xp_exc are identical for strawberries and ranges from 0.41 to 0.51. This approach allows the calculation of bulk dried tissue δ¹⁸O values without the knowledge of cellulose δ¹⁸O 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 δ¹⁸O values are assessed by an offset (ε_cp) to the cellulose δ¹⁸O values.

Model parameter selection

To find the best values of the key model parameters for the prediction of strawberry bulk dried tissue δ¹⁸O 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) f_xylem and p_xp_ex values reported in literature for leaf water and cellulose from various species that were averaged, (ii) values averaged for leaves (f_xylem) and berries (p_xp_exc) of berry producing plants, or (iii) values for leaves (f_xylem) and berries (p_xp_exc) 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 δ¹⁸O for different species and averaged these values (0.22 for f_xylem^{36,37,38,39,40} and 0.40 for p_xp_ex^{45,46,47,48,49}) (Table 1). For berries (average of the values of raspberries and strawberries) the mean leaf-derived f_xylem value was 0.26 and the value for p_xp_exc was determined to be 0.46 (Table 1) (data derived from Cueni et al. in review). For strawberry plants, the leaf-derived f_xylem value we used was 0.30, and the value for bulk dried tissue (p_xp_exc) was determined to be 0.46 (Table 1) (data derived from Cueni et al. in review). Since the p_xp_exc values of different berry species did not differ, this resulted in a total of six different model input parameter combinations.

Environmental model input data selection

In 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 δ¹⁸O 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 season^46,53, and precipitation δ¹⁸O values of rain-events prior and during the growing season^51,54,55 have been shown to shape plant tissue water and organic compound δ¹⁸O 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).

Two precipitation isotope data products were compared (Table 2): (1) The mean monthly precipitation δ¹⁸O 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)³². This source provides values for individual months and years based on station coordinates³². Both data sets were on the one hand used for the precipitation δ¹⁸O input data of the model, and also to extrapolate the vapor δ¹⁸O 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)²⁹ and the E-OBS gridded dataset by the European Climate Assessment & Dataset (Version 22.0e)³⁰ (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 date^46,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 averages⁷⁶, 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 δ¹⁸O 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 samples

Using the plant physiological model described above, we calculated the strawberry bulk dried tissue δ¹⁸O 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 (f_xylem and p_xp_ex/p_xp_exc, Table 1), yielding model results to be evaluated against the measured reference samples. Our approach can be described with the following equation:

where δ¹⁸O plant is the simulated δ¹⁸O 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 f_xylem and p_xp_ex/p_xp_exc we on the one hand used the values proposed for leaves by literature (for p_xp_ex), 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 p_xp_exc). 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 model⁷⁷, 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 analyses

Statistical analyses were done using the statistical package R version 3.5.3⁷⁹. The relationships of the range δ¹⁸O 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 δ¹⁸O bulk dried tissue values of the authentic reference samples (n = 154) by calculation of the root mean squared error (RMSE).

Calculation of prediction maps

Prediction 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 δ¹⁸O 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 δ¹⁸O values of berries grown in July 2017. For this we used the average berry model input parameters (f_xylem and p_xp_exc, 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 δ¹⁸O values from OIPC as an average from April, May and June, and vapor δ¹⁸O values calculated from OIPC precipitation δ¹⁸O 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 δ¹⁸O 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 δ¹⁸O value and the predicted map δ¹⁸O value for each pixel of the map. The places (pixels) that are predicted to have the same δ¹⁸O 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 δ¹⁸O value of the suspected sample could be assessed. This means that the bigger the difference between the simulated δ¹⁸O value and the sample value, the lower the probability of provenance of the sample. In other words, a difference between the sample’s δ¹⁸O value and the predicted δ¹⁸O 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).

Source: Ecology - nature.com