key: cord-0028675-glwabfsf authors: Kramer, R. T.; Kinaston, R. L.; Holder, P. W.; Armstrong, K. F.; King, C. L.; Sipple, W. D. K.; Martin, A. P.; Pradel, G.; Turnbull, R. E.; Rogers, K. M.; Reid, M.; Barr, D.; Wijenayake, K. G.; Buckley, H. R.; Stirling, C. H.; Bataille, C. P. title: A bioavailable strontium ((87)Sr/(86)Sr) isoscape for Aotearoa New Zealand: Implications for food forensics and biosecurity date: 2022-03-16 journal: PLoS One DOI: 10.1371/journal.pone.0264458 sha: 4ad555f743c3b7bec6426bf6afb6d5105b1d1f3b doc_id: 28675 cord_uid: glwabfsf As people, animals and materials are transported across increasingly large distances in a globalized world, threats to our biosecurity and food security are rising. Aotearoa New Zealand is an island nation with many endemic species, a strong local agricultural industry, and a need to protect these from pest threats, as well as the economy from fraudulent commodities. Mitigation of such threats is much more effective if their origins and pathways for entry are understood. We propose that this may be addressed in Aotearoa using strontium isotope analysis of both pests and products. Bioavailable radiogenic isotopes of strontium are ubiquitous markers of provenance that are increasingly used to trace the origin of animals and plants as well as products, but currently a baseline map across Aotearoa is lacking, preventing use of this technique. Here, we have improved an existing methodology to develop a regional bioavailable strontium isoscape using the best available geospatial datasets for Aotearoa. The isoscape explains 53% of the variation (R(2) = 0.53 and RMSE = 0.00098) across the region, for which the primary drivers are the underlying geology, soil pH, and aerosol deposition (dust and sea salt). We tested the potential of this model to determine the origin of cow milk produced across Aotearoa. Predictions for cow milk (n = 33) highlighted all potential origin locations that share similar (87)Sr/(86)Sr values, with the closest predictions averaging 7.05 km away from their true place of origin. These results demonstrate that this bioavailable strontium isoscape is effective for tracing locally produced agricultural products in Aotearoa. Accordingly, it could be used to certify the origin of Aotearoa’s products, while also helping to determine if new pest detections were of locally breeding populations or not, or to raise awareness of imported illegal agricultural products. Introduction not be distinct enough. This issue is most important in Aotearoa and the greater Pacific where islands and small land masses do not experience the same continentality effects that lead to distinct isotopic variation in other regions of the world. The other light isotopes, δ 13 C and δ 15 N, are primarily used to look at variations in diet for humans and animals and farming practices [22, [37] [38] [39] [40] [41] [42] . Singularly, these isotope tracers have limited application to geographical origin prediction but combining them with other isotopes can assist with provenancing efforts. This is because δ 13 C and δ 15 N vary predictably with the climate, dietary resources, and long-term land-use effects as they differentiate between C 3 and C 4 plants, protein consumption in the diet, and identifying the use of conventional and organic fertilizers in farming systems [37, 40, 43] . Lastly, all light isotopes are susceptible to variation introduced through the global supermarket, where nonlocal products may be incorporated into modern human diets and blur the geochemical signature within the sample [19, 37, 44, 45] . In addition to the light stable isotopes used for provenancing, trace elements (chemical elements found in low concentrations of less than 100 parts per million (ppm) in the mineral matrix) are also not used to predict origins but can be used for samples of known origin to determine if they have similar or distinguishable chemical profiles from other, potentially fraudulent, samples [46, 47] . Chemical profile comparisons can be based on just a few key elements or on a whole suite of elements depending on the sample and reference materials available [47] . For example, trace elements could identify fraudulent food products, like wine or tea, by comparing the chemical fingerprint of the suspicious sample to the authentic product [20, 21] . The caveat is that this requires a priori knowledge about the chemical profile for all potential regions of origin to determine if the sample of interest classifies into a particular region or group [20] . While feasible, this would require an extensive database to store the chemical fingerprints for every type of reference material that may need to be provenanced. Furthermore, Koffman et al. [48] found that trace elements from Aotearoa sediments lack the regional variability observed when using lead, neodymium, and strontium isotopes, but trace elements have successfully been used to demonstrate regional variability in Aotearoa soils at multiple scales [49, 50] . The analysis of the radiogenic strontium isotope ratio, 87 Sr/ 86 Sr, is an alternative approach that has the potential to be broadly applicable without the need for situation-specific research and development. Internationally, this isotope system has been a key investigative method used to predict the region-of-origin for plants, animals, insects, and other biological materials [18, [51] [52] [53] [54] [55] [56] . Ecosphere 87 Sr/ 86 Sr variation primarily reflects the underlying geology but also other processes including chemical weathering, alluvial and fluvial erosion, soil processes, and aerosol deposition (sea salt, volcanic ash, dust, loess) [51, 55, 57] . Strontium in the geosphere enters the tissues of biological organisms through their uptake of water and ingestion of dietary resources. This integrated 87 Sr/ 86 Sr fraction is referred to as biologically available or "bioavailable". Any minor isotopic fractionation that occurs in the environment is corrected accordingly during analysis [51, 52, 55, 58, 59] . Therefore, biological tissues usually retain the "isotopic fingerprint" of the local ecosystem from which dietary resources were obtained [58] . Different biological tissues (plant leaves, hair, nails, teeth, and bone) form and remodel at different rates and the integrated 87 Sr/ 86 Sr values in these tissues should reflect where organisms lived at different time frames of their lives [60] [61] [62] . Strontium isotope analysis as a tool in forensic tracing has received limited attention to date in Aotearoa mostly due to the cost of analysis when scaled to very large numbers of samples, and the lack of a bioavailable 87 Sr/ 86 Sr baseline. Recent efforts have gone some way to tackling these issues, with technical advances to address the former [63] and with Duxfield et al. [64] aggregating published geological data to create a 87 Sr/ 86 Sr baseline. The recent baseline [64] summarized the variation in geological 87 Sr/ 86 Sr for Aotearoa and did not predict the bioavailable 87 Sr/ 86 Sr produced through the complex environmental system that 87 Sr/ 86 Sr cycles through. Duxfield et al. [64] used a case study to compare plant bioavailable 87 Sr/ 86 Sr values to their geologic baseline and found that the values fell within the expected ranges based on the underlying lithology but demonstrated smaller 87 Sr/ 86 Sr value ranges than the geological unit they were associated with. These results are most likely due to the geological 87 Sr/ 86 Sr baseline not considering the exogenous sources of 87 Sr/ 86 Sr in the biosphere that are influenced by climatic, atmospheric, and environmental conditions. Also, some of the published geological data they used are based on 87 Sr/ 86 Sr compositions for individual mineral phases rather than the bulk rock, which homogenizes mineral signatures. Similar to baselines, isoscapes are models that predict the spatial distribution of isotopes by incorporating various sources that contribute to the bioavailable isotope "pool" of a region [65] [66] [67] [68] [69] [70] [71] [72] . Strontium isoscapes have been produced at the continental scale for North America [73] , Central America [74] , Africa [75] , and Europe [76] , at regional scales for Australia [77] , the Caribbean [78] , China [79] , the Netherlands [80] , France [81] , southwest Sweden [82] , the Modena [83] , and South Korea [84] , as well as for island locales such as mainland UK [85] and Ireland [86] . For strontium isotopes, the leading approach to create regional and global isoscapes [51,56,87-89] uses bioavailable 87 Sr/ 86 Sr values of georeferenced samples (soil, plants, and small local organisms) and a machine-learning framework that considers contributions from the ecosphere to predict bioavailable 87 Sr/ 86 Sr values across the region of interest. In this paper, we introduce the first bioavailable 87 Sr/ 86 Sr isoscape model for Aotearoa constructed using a machine-learning approach [51, 88] and the best available dataset of geospatial predictors. Specifically, the isoscape construction uses a random forest (RF) model framework that considers a variety of climatic, atmospheric, and environmental variables that may contribute to the spatial distribution of bioavailable 87 90, 91] . To compensate for this, the sampling strategy of this study targeted plants of varying root depths (shallow, medium, and deep) at each sampling location to ensure that any intra-site 87 Sr/ 86 Sr variability was captured. Shallow roots are defined as slender, branched, fibrous or creeping roots that grow close to the surface. Medium roots include larger plants whose roots penetrate one to two meters below the surface soil. Deep-rooted plants have tap roots that consist of a primary root that penetrates deep into the soil at depths greater than three meters. Using the RF-based methodology allowed us to explain the main environmental influences on the 87 Sr/ 86 Sr isoscape and we then validated the use of for provenancing using cow milk samples collected from farms across Aotearoa. Plant sample collection sites were chosen based on a strict set of inclusion criteria to circumvent potential anthropogenic and natural Sr contaminants. These criteria avoided farms, pastures, drainage ditches, and other human-made earthworks, and being restricted to public lands with road-access. Riparian areas, land in flood areas and along watercourses, were avoided to avoid any mixing of exogenous Sr sources [92] . Also, the topographic relief and elevation of the regions limited where and when samples could be taken. These factors greatly limited the distribution of the samples (Fig 1) . Potential sample stops were planned from an ArcGIS shapefile layer with a highway road feature class overlaid on a geological map to mark where there was a geological change within each area along the route. Leaves were collected from plants of each root-type (shallow, medium, and deep) located within a 5 m diameter of one another at each sample site (when available) to capture a higher resolution of 87 Sr/ 86 Sr variation for the single locality. Samples were stored in labeled, plastic sample bags and dehydrated later the same day in a clean SunBeam food dehydrator before being stored in new clean and labeled, plastic sample bags. Plant type, root-depths, and other metadata are provided in the Supplementary Information (S3 File). Topsoil samples (0-20 cm) were collected by GNS Science in the Nelson, Otago, and Southland regions (Fig 1) . Samples were collected by hand auger and the sub-2-mm portion was retained after sieving and drying at 40˚C [49, 93] . Additional soil and plant 87 Sr/ 86 Sr data were provided by PWH and KFA for the Christchurch, Bay of Plenty, Auckland, and Northland regions. To date, very few modern local human and animal bioavailable 87 Sr/ 86 Sr data have been published for Aotearoa and only eight (two sheep and six humans) were available to include in this study (Fig 1) . Multiple analytical methods were utilized to analyze the plant and soil samples and generate the 87 Sr/ 86 Sr data. GNS topsoil (n = 71) and collected plant samples (n = 185) were prepared under the supervision of the lab technician in the Centre for Trace Element Analysis, University of Otago, Dunedin. GNS topsoil subsamples of 1g were leached in 2.5 mL of 1 M ammonium nitrate (NH 4 NO 3 ) solution and agitated overnight to extract the bioavailable strontium fraction [53, 81, 92, 94] . Then, the topsoil samples underwent Sr separation column procedures as described previously [76, 95, 96] , while the collected plant samples utilized an automated ion-exchange chromatography method using a 3-ml column (ESI, Part number CF-MC-SrCa-3000) filled with DGA (diglycolamide) resin (TrisKem International, Bruz, France) as described by Wijenayake [97] . Strontium isotope measurement of the collected plant and GNS topsoil samples was conducted at the Centre for Trace Element Analysis, University of Otago using a Nu Plasma-HR MC ICP-MS (Nu Instruments Ltd., UK) following previously reported protocols [76, [95] [96] [97] [98] [99] . Repeated measurement of the NIST SRM 987 strontium isotope reference material (sourced from the National Institute of Standards and Technology (NIST), USA) and the HPS inhouse strontium isotope standard, that bracketed every six samples, were used to monitor the accuracy and external precision of the measurements. We obtained average values of 0.71025 ± 0.00002 (2 SD, n = 70) for NIST SRM 987 and 0.70761 ± 0.00009 (2 SD, n = 58) for HPS in very good agreement with previously reported values. Any instrumental mass fractionation present was corrected for using the exponential mass fractionation lay by normalization to 86 Sr/ 88 Sr = 0.1194. Procedural blanks were run with each batch of 6 samples and all yielded negligible Sr levels of < 250 pg. The additional soil and plant samples (n = 126) utilized Sr separation column procedures as described by Pin & Bassin [100] . Sr isotope measurement of the soil samples from [101] was conducted using a Nu Plasma-HR MC ICP-MS (Nu Instruments Ltd., UK) at the Centre for Trace Element Analysis, University of Otago, Dunedin (standard = NIST SRM 987, n = 19, 87 Sr/ 86 Sr 0.710274 ±0.000023 (2 SD)); and the additional soil and plant samples [102] measured using Isotopx Phoenix TIMS (thermal ionization mass spectrometry) at the University of Adelaide (NIST SRM 987, n = 7, 87 Sr/ 86 Sr = 0.710245 ± 0.000008 (2 SD). Further details regarding preparation and analysis for all plant and soil samples are available in the Supplementary Information (S1.1-5 in S1 File). Altogether, 414 bioavailable 87 Sr/ 86 Sr data (Fig 1) were used to construct the strontium isoscape, specifically comprising data from 314 plants (182 unique locations), 92 topsoils (84 unique locations), and eight mammals (two unique locations). This study follows the methodology established by Bataille and colleagues [51,93] using a RF model framework. Model construction used bioavailable, geo-referenced 87 Sr/ 86 Sr data gathered for Aotearoa, many geomatic auxiliary variables in the form of global rasters (gridded matrix of cells) used as covariates, and an existing process-based bedrock 87 Sr/ 86 Sr model [87, 88] to produce the final predicted 87 Sr/ 86 Sr isoscape model. The RF R-script provided by Bataille et al. [51, 88] optimizes the regression model with the root mean square error (RMSE) and uses five 10-fold repeated cross-validations. Additionally, the R-script includes using the Variable Selection Under Random Forest (VSURF) package [103] to identify relevant and highly predictive variables [51, [87] [88] [89] . The relationships between the variables selected by the VSURF function and the bioavailable 87 Sr/ 86 Sr variability were assessed using variable importance purity measure and partial dependence plots [51] . The auxiliary variables were obtained from various sources summarized in S2 Table 1 in S2 File and represent geological, climatic, and environmental variables that may influence bioavailable 87 Sr/ 86 Sr variability. Most variables were available as global rasters that were trimmed to an Aotearoa extent and projected to the New Zealand Transverse Mercator 2000 coordinate system. To assess the uncertainty of the isoscape, we used quantile RF regression model to generate quartile-1 and quartile-3 models using the log-transformed 87 Sr/ 86 Sr values. We use a logtransformation because the 87 Sr/ 86 Sr data have a positively skewed distribution. Then, the quartile-1 model was subtracted from the quartile-3 model (i.e., Q3-Q1) to generate a final interquartile range raster at a resolution of 1 km. Dairy farms throughout the country that contributed cow milk samples agreed to a controlled feeding regime where all cattle were expected to be pasture-fed on-site and were not provided with supplementary feed options [97] . Therefore, the assumption was that the 87 Sr/ 86 Sr values of the cow milk would reflect the underlying geology and local atmospheric conditions with no contamination from exogenous 87 Sr/ 86 Sr sources. Details regarding the isotopic preparation and analytical methods for the cow milk samples are provided in the Supplementary Information (S1.5 in S1 File). To assess the validity of the bioavailable 87 Sr/ 86 Sr isoscape developed, we predict the regionof-origin for cow milk samples and then assess the accuracy of that assignment relative to their actual origin. Firstly, origin predictions are computed using assignR [104] , which operates in a semi-parametric Bayesian framework to calculate the posterior probability of the sample belonging to each cell within the isoscape raster. Statistically, assignR uses a maximum likelihood assignment model and employs Bayes Theorem to calculate the probability that a sample originates from a geographic location given the isotopic signature of the sample [105] . The statistical methodology of the assignR code [104] is adapted from Wunder [71] and Vander Zanden et al. [106] . To facilitate visualization, we created maps displaying the top 33% by area of the predicted surface for each sample (SI4) with the top 33% areas coded as 1 and other areas coded as 0. Once the top 33% probability maps are produced, they are brought into ArcGIS Pro. Then, we calculated the accuracy of the nearest predicted geographic origin (part of the top 33% probability) for each sample in terms of how close they were in km to the actual place of origin for the cow, using the Measure Distance tool in ArcGIS Pro. We also measured the distance from the nearest predicted cell for the top 20% and top 10% probability surfaces to assess the trade-offs between accuracy and precision for these arbitrary thresholds. We do not provide the maps for the top 20% and 10% predictions, but they can be provided upon request. The bioavailable 87 Sr/ 86 Sr data for the plant and soil samples are available in the Supplementary Information (S3 File). The combined values for plants, soils, and local mammals demonstrate comparable distributions, with Quartile 1 = 0.70738 ± 0.00002, median = 0.70832 ± 0.00003, and Quartile 3 = 0.70897 ± 0.00108 (Fig 2) . The animal substrate samples display a tighter range of 87 Sr/ 86 Sr values compared to the plants and topsoils, with soil displaying the largest variability in 87 Sr/ 86 Sr, but 90% of all data fall within values of 0.70500 and 0.71250 (Fig 2) . The final 87 Sr/ 86 Sr isoscape for Aotearoa uses the best performing RF model that considered all auxiliary variables (S2 Table 1 in S2 File) and 414 bioavailable 87 Sr/ 86 Sr values (S3 File). The RF regression model produces a map demonstrating the mean 87 Sr/ 86 Sr prediction (R 2 = 0.53, RMSE = 0.00098) ranging from 0.70567 to 0.71118, for the entire country including the Chatham Islands (Fig 3A) . The accompanying interquartile range raster ( Fig 3B) demonstrates a standard error ranging from 0.0001 to 0.002 for the country. The VSURF package selected 11 predictive variables ( Fig 4A) . These include dust and sea salt aerosol deposition (r.dust, r.ssa, r.ssaw), geological attributes (r.age, r.toprock, r.GNSagemax), temperature (r.mat), elevation (r.elevation), and soil characteristics (r.ph, r.clay, r.nitrogen). Initially, the VSURF selected 10 features, but when r.nitrogen was forcibly added to the RF function, the amount of variation explained increased by~2%. Therefore, r.nitrogen was included in the final model regression. The dominant predictive variables are the rock type classification (r.toprock) and soil pH in H 2 O solution (r.ph) based on their Percent Increased Mean Squared Error (%IncMSE) values ( Fig 4A) . The higher the %IncMSE value, the more important the predictor variable [107] . The n-fold cross validation explains 53% of the variation in the bioavailable 87 Sr/ 86 Sr model, with an RMSE of 0.00098 over the dataset (Fig 4B) . The uncertainty appears uniform across the prediction range ( Fig 4B) . Partial dependence plots are used to examine the association between the predictors and bioavailable 87 Sr/ 86 Sr (Fig 5) . These illustrate that 87 Sr/ 86 Sr values increase with increasing soil pH (r.ph) and elevation (r.elevation). On the other hand, inverse relationships are apparent where 87 Sr/ 86 Sr values increase with decreasing mean annual temperature (r.mat), principal surface lithology type (r.toprock), rock age (r.age, r.GNSagemax), and clay content (r.clay). The "toprock" raster ( Fig 5) consists of 67 different lithological classifications (detailed in S2 Table 2 in S2 File) where the first few categories have the highest 87 Sr/ 86 Sr values: 1 = floodplain alluvium; 2 = volcanic ashes older than Taupō pumice; 3 = Loess; 4 = Sandstone; Units of measurement are provided on the x-axis for each variable. For r.toprock, the x-axis represents the numbered categories of rock type (S2 Table 2 [73, 87, 88] . The remaining predictive variables, dust and sea salt aerosol deposition (r.dust, r.ssa, r.ssaw) and nitrogen content (r.nitrogen), display non-linear relationships with the 87 Sr/ 86 Sr values. Generally, we found that cow milk with 87 Sr/ 86 Sr values near 0.70800 were difficult to assign because most of the isoscape values across Aotearoa fall in a narrow band between the range of 0.70750 and 0.70930. When a sample is within this interval, the probability surfaces show broad regions with high probabilities of origin corresponding to low assignment precision. Conversely, the model performs well when predicting origin for samples that are below 0.70750 or above 0.70930 because these values are much less common on the Aotearoa isoscape allowing for better precision. We produced prediction maps showing the top 33% of isoscape raster cells with the highest posterior probability for each cow milk sample (n = 33) available in the Supplementary Information (S4 File). The cow milk 87 Sr/ 86 Sr and the predicted 87 Sr/ 86 Sr from the isoscape show a good correlation (R 2 = 0.52) when compared with one another (Fig 6) . When error bars are included, all cow samples fall within the 95% confidence interval of the linear regression line on the actual versus predicted bioavailable 87 Sr/ 86 Sr value plot (Fig 6) , except sample Cow 7. The maximum likelihood assignment model was quite accurate and all but one milk sample had high probability cells (top 33%) located within proximity (average of 7.05 kilometers) to (Tables 1 and 2 ). We also measured the distance between the nearest probability cell and the dairy farm of actual origin for the top 20% and top 10% probability surfaces to assess the trade-off between precision and accuracy when different probability thresholds are used (Tables 1 and 2 ). While the accuracy of most origin predictions is high using the top 33% probability threshold, the precision is lacking for most samples. Conversely, the top 20% and top 10% probability thresholds produce precise origin predictions (i.e., they predict fewer regions of potential origin), but the accuracy (distance to known origin) decreases (Tables 1 and 2) . Cow milk samples originating from the Canterbury region produce the most accurate predictions (Tables 1 and 2 ). For all Canterbury cow milk samples, the prediction maps have high probability cells within 0 to 1.56 km of the actual location. The Waikato region milk samples Understanding the geology and climate allows improved modeling of Sr isotope cycles from the geosphere into the biosphere for provenancing materials. Accordingly, the covariates used in the RF model here account for geological and atmospheric conditions in Aotearoa (S2 Table 1 Fig 3) across Aotearoa. These values reflect the variable composition of bedrock for the country, which has a rich geological history [108] that has resulted in areas with distinct mineralogy and rock ages, both factors which strongly affect strontium and consequently bioavailable 87 Sr/ 86 Sr [109] . Aotearoa has a compositionally diverse basement geology with formation ages between c. 500-100 million years ago (Ma) (Fig 7) . The modern-day distribution of geology reflects terrane accretion and c. 23 Ma of tectonism along the strike-slip system, the Alpine Fault, which bisects the country [108, 110, 111] . The Alpine Fault juxtaposes terranes which were once contiguous for hundreds of kilometers, resulting in rocks of terrane affinity being present at either end of the country. belonging to the Paleozoic Buller and Takaka terranes of the Indo-Australian plate. The rocks of these terranes are broadly comprised of metasediments and intruding plutonic igneous rocks, called batholiths, that formed along the previously active margin of Gondwana between 500-100 Ma [112, 113] . Many of these rocks are of igneous origin or derivation, and therefore tend to have very low 87 Sr/ 86 Sr values despite their older age [114, 115] , though some plutonic igneous rocks are known to have higher 87 Sr/ 86 Sr values owing to their derivation from melting ancient crust [116] . The lowest predicted bioavailable 87 Sr/ 86 Sr isoscape values in the North Island also partially correspond with basement geology-some areas of low predicted 87 Sr/ 86 Sr isoscape values correspond with the Dun Mountain-Maitai and the Murihiku terranes (Fig 7) . The Dun Mountain-Maitai Terrane is composed of ultramafic rocks overlaid by plutonic and volcanic sequences, that are covered by sedimentary deposits, while the Murihiku Terrane consists largely of sedimentary and volcanic materials [118] . However, in the central North Island bioavailable strontium ratios also appear to be heavily influenced by Quaternary volcanic deposits relating to eruptives from the Taupō Volcanic Zone (S2 Fig 1 in S2 File) . These eruptive rocks include ignimbrites and volcanic ashes that date volcanism from around 1.6 Ma to present day [119] . Therefore, these deposits and the ashes that predate them are young and are expected to display low 87 Sr/ 86 Sr values. The area immediately surrounding Mount Taranaki in the westernmost region of the central North Island (Fig 7) has been active since c. 130 thousand years ago (Ka) [120] and is also comprised of volcanic rocks that display the low 87 Sr/ 86 Sr values characteristic of younger igneous formations [55] . The partial dependence plots (Fig 5) demonstrate this relationship between 87 Sr/ 86 Sr values and age as well. Conversely, the highest 87 Sr/ 86 Sr values (0.70930-0.71120) of Aotearoa are predicted for areas comprised of the Mesozoic Caples, Rakaia, and Pahau terranes which stretch from the lower South Island along the eastern portion of the North Island (Fig 7) . These terranes have mudstone/sandstone protoliths that are variably metamorphosed, but generally increase in metamorphic grade the closer they are to the Alpine fault [121, 122] . The Caples Terrane consists of volcaniclastic feldspathic metasedimentary rocks [123] . The Rakaia and Pahau terranes contain variably metamorphosed quartzofeldspathic sedimentary rocks [124] . These rocks are overlain by aerially expansive alluvium and loess deposits in the central east portion of the South Island, around the Canterbury Plains (Fig 7 and S2 Fig 1 in S2 File) . The elevated 87 Sr/ 86 Sr values for these regions (Fig 7) most likely result from the more felsic (rich in silica) content of the rocks [64] . The juxtaposition of differing geologies along the Alpine Fault also has clear effects on the distribution of bioavailable 87 Sr/ 86 Sr values in Aotearoa. This fault is primarily a strike-slip system, however uplift on the southeastern side of the fault has resulted in the formation of the Southern Alps. This mountain chain runs along the Alpine Fault and marks the distinction between high and intermediate 87 Sr/ 86 Sr values in the model (Fig 7) . The elevation partial dependence plot (Fig 5) shows that 87 Sr/ 86 Sr values increase with increasing elevation. In addition, erosional forces associated with uplifted areas have generated massive amounts of gravel and dust that are older and have high 87 86 Sr, atmospheric variables clearly interact with geology to produce a complex isoscape. Despite its small size, Aotearoa's climate varies vastly with a warm subtropical climate in the Northland region, a cool temperate climate throughout most of the country, and alpine weather in the mountainous West Coast and Southland regions of the South Island [125] . The mean annual temperature plot (Fig 5) demonstrates that temperature has a significant effect on bioavailable 87 Sr/ 86 Sr in Aotearoa where higher 87 Sr/ 86 Sr values occur in areas with lower temperatures. This pattern correlates with geological and latitudinal differences between the North and South Islands. The former experiences warmer weather and has more volcanic and low-Sr lithologies (Fig 7) , while the latter experiences colder weather and the highest 87 Sr/ 86 Sr values occur along the Southern Alps and the Canterbury plains. Precipitation is another atmospheric variable that plays a key role in bioavailable 87 Sr/ 86 Sr distributions. The Southern Alps experience some of the highest recordings for mean annual precipitation in Aotearoa, measuring up to 6,000 mm/year, while the land to their east records the lowest rainfall at 250-1,000 mm/year [125, 126] . Increased precipitation on this coastal setting probably facilitates a high rate of marine aerosol deposition that is characterized by lower 87 Sr/ 86 Sr signatures, as well as historically high deposition of volcanic ash from surrounding volcanic centers particularly in the northwestern South Island and the North Island [55] . This is reflected in the partial dependence plot for the sea salt aerosol deposition where higher rates of deposition are associated with lower 87 Sr/ 86 Sr values (Fig 5) . However, the lack of bioavailable 87 Sr/ 86 Sr samples from these regions does limit the accuracy of the isoscape here (Fig 7) . No cow milk samples were obtained from these regions. Using the first bioavailable 87 Sr/ 86 Sr isoscape developed here for Aotearoa, this study sought to validate the use of strontium isoscapes for provenancing research, such as food product verification, by predicting the origin for cow milk sampled from around the country. Fig 6 illus trates that all cow milk samples, with one exception (Cow 7), fall within the 95% confidence interval of the line of best fit. This indicates that their measured 87 Sr/ 86 Sr values reflect the predicted values obtained from the bioavailable 87 Sr/ 86 Sr isoscape, illustrating that the bioavailable 87 Sr/ 86 Sr isoscape performs well and has the potential to be utilized in provenancing applications. The relatively high regression residual for sample Cow 7 might be due to the addition of non-local feed to the cow's diet including health supplements, supplementary feed from other Aotearoa regions including ryegrasses (with high calcium content), imported palm kernel expeller (PKE), or even transportation of the cow between regions in Aotearoa that could introduce foreign 87 Sr/ 86 Sr values. However, dairy farms that contributed cow milk samples [97] had enforced a controlled feeding regime where all cattle were expected to be pasture-fed on-site and were not provided with supplementary feed options. In Aotearoa, livestock and raw cow's milk are commonly transported inter-regionally. Prior to going through the pasteurization and production processes, and depending on processing plant availability or maintenance, milk may be transported some distance, between the west and east coasts, particularly in the South Island. The outlier sample, Cow 7, may therefore represent a cow that 1) was not kept on a strict diet or 2) had been transported from another region in Aotearoa. The accuracy and precision of predictions was measured as the distance from their known place of origin to the nearest predicted area of potential origin using different probability quantile thresholds (10%, 20%, and 33%). As expected, the top 33% threshold produced the most accurate (i.e., closest to known origin) predictions, but those predictions covered more potential areas and were less precise (Tables 1 and 2 and S4 File). On the other hand, the top 20% and top 10% thresholds provided more constrained potential region-of-origin predictions, making them more precise, but their accuracy decreased (Tables 1 and 2 ). Specifically, the average distance away from the place of known origin was 7.05 km for the top 33% threshold, 16.99 km for the top 20%, and 30.73 km for the top 10% threshold (Table 2 ). This showed that there was a trade-off between precision and accuracy for prediction outputs. In the Supplementary Information (S4 File), we provide prediction maps of the top 33% probability quantile for each cow milk sample. In total, 73% (24 out of 33) of milk samples fell into the intermediate 87 Sr/ 86 Sr range with values between 0.70750 and 0.70930 (Figs 6 and 7) . Milk that fell into this range originated from the Waikato, Canterbury, Southland, Northland, Manawatu-Wanganui, and Nelson regions (Figs 1 and 6) . As indicated by the isoscape, there is high redundancy for bioavailable 87 Sr/ 86 Sr across Aotearoa. When samples had 87 Fig 6, and S4 File) . When 87 Sr/ 86 Sr values did not fall into the low or high ranges, the accuracy of origin predictions PLOS ONE decreased, meaning that the distance away from the known place of origin tended to be larger because the assignment model predicted many potential regions of origin with similar bioavailable 87 Sr/ 86 Sr values. In these cases, utilizing a second isotope system, such as δ 2 H and δ 18 O, may help to constrain the predicted region-of-origin [31, 37] . In terms of accuracy, Canterbury region cow milk samples (n = 8) produced the most accurate origin predictions across all probability thresholds with the smallest average distance between the known origin and the nearest predicted cell ( Table 2 and S4 File). This strongly implies that the Canterbury cows have a local diet and do not have foreign feed introducing exogenous 87 Sr/ 86 Sr values. Conversely, Southland region milk samples (n = 11) produced the least accurate origin predictions across all thresholds (Tables 1 and 2 ). Fig 6 shows that the majority of the Southland cow milk samples have observed values below the 1:1 black line, implying that their actual 87 Sr/ 86 Sr values are higher than the predicted 87 Sr/ 86 Sr isoscape values. As mentioned previously, much of the South Island was glaciated about~26,000 to~18,000 BP and glacially fed rivers, specifically the Clutha and Mataura rivers and their tributaries, carried sediment from the Wakatipu Valley in the far western portion of the Otago region into the Southland and eastern Otago regions [48] . The fluvial transportation of sediments from the Wakatipu Valley may be a major contributor to the elevated 87 Sr/ 86 Sr values of the Southland cow milk samples. Furthermore, the low accuracy may result from the increased geologic diversity in Southland versus Canterbury, where the former is comprised of several terranes, ranging from mafic to felsic, and the former is dominated by the Rakaia Terrane (Fig 7) . Even though the Southland cow milk samples fell below the 1:1 line (Fig 6) , many of the samples (Cow 19-23, 26, and 27) produced accurate predictions less than 5 km away from their known origin and three other Southland samples (Cow 24, 25, and 28) were only 13 to 19 km off their targets using the top 33% threshold (Table 1) . One sample, Cow 29, had the highest 87 Sr/ 86 Sr value of the Southland group and was the least accurate where the nearest predicted cell was 86-123 km away from the cow's known place of origin (Table 1 ). Since the other milk samples from the Southland region produce accurate to moderately accurate predictions, we do not suspect that the elevated 87 Sr/ 86 Sr value of sample Cow 29 was due to its proximity to the coast (i.e., sea salt or dust aerosol deposition) or the introduction of foreign feed. Instead, it is likely that Cow 29 had been transported down to the Southland region from the Canterbury region which is more consistent with its elevated 87 Sr/ 86 Sr values. Milk samples were obtained through farms on a voluntary basis, so determining whether this sample had actually been transported from Canterbury was not possible. Apart from food science applications, the bioavailable 87 Sr/ 86 Sr isoscape can assist with identifying the origin of new pests, unidentified materials, illegal agricultural products, and origin mislabeling. Specifically, 87 Sr/ 86 Sr can help biosecurity organizations within Aotearoa determine if a pest encountered post-border represents a newly introduced pest or an established population bearing local 87 Sr/ 86 Sr values [1, 99, 100, 127, 128] . For the key pastoral industry, Ferguson et al. [127] estimated that invasive invertebrates alone cost the sector between $1.7 to $2.3 billion NZD annually. New pests and diseases are introduced to Aotearoa through shipping and air traffic pathways [128] [129] [130] [131] . One of the highest risk pests for Aotearoa is the exotic brown marmorated stink bug (Halymorpha halys) because it feeds on over 300 species of plants and infestation can ruin entire crops [7] . Recent attempts to define the origins of this invasive stink bug used δ 2 H and δ 18 O isotopes from their wings to determine whether a bug detected post-border represented a recently introduced foreign bug or an established population [2] . Holder et al. [2] concluded that the distinct δ 2 H and δ 18 O isotope signature of their wings suggested a cooler climate origin, supporting evidence that the specimen was not from a locally breeding, southern hemisphere summer population. Currently, there are no rapid response tools to assess the provenance of these pests when they are introduced which limits the ability to develop effective measures to prevent their arrival in Aotearoa. This leaves agencies poorly informed as to the actual risk, with difficult decisions about embarking on expensive responses that may otherwise have to assume the presence of a locally established population. Recently, Murphy et al. [63] developed a new mass spectrometric method to analyze 87 Sr/ 86 Sr in low-Sr samples, such as a single insect, in less time than traditional methods. This new analytical method enables this technology to be used within the limits for biosecurity, commanding very short time frames and availability of very little biological material. Now, coupled with this 87 Sr/ 86 Sr isoscape, well supported decisions based on probabilities can be made as to the likely Aotearoa or offshore origin of foreign pests, which is key to determine the most efficient response to prevent their establishment [2] . Prior to this point this had not been possible, anywhere. The identification and removal of foreign pests protects crops from foreign diseases and safeguards the future of endemic flora and fauna species throughout the country. This 87 Sr/ 86 Sr-focused approach can be applied globally to determine the local or foreign nature of pests on borders and ports of entry. However, this would require that areas of interest construct bioavailable 87 Sr/ 86 Sr isoscapes to compare their samples against. The Aotearoa 87 Sr/ 86 Sr isoscape produced by this study could be used alongside other lines of evidence (e.g., isotopic systems, chemical fingerprints), and strontium baseline models from other regions of the world to test whether "Made in New Zealand" food products display expected 87 Sr/ 86 Sr signatures or if they are fraudulent products. Furthermore, law enforcement agencies could use the isoscape to predict the region of 'growing areas' or determine country of origin using bioavailable 87 Sr/ 86 Sr from confiscated illicit drugs, like marijuana or heroin [132, 133] . Additionally, although Aotearoa does not encounter many unidentified human remains, if they were recovered and traditional identification methods failed to produce a positive identification, 87 Sr/ 86 Sr values from the unidentified person's teeth could help predict their region or country of origin [19, 25, 26, 37, 65, 74, 77, 134, 135] . The method of provenancing used in this study is useful because it provides a visual representation of the predicted probability surface that allows for ease of communication and comprehension. Though the predicted maps need to be generated by a specialist, once created, a non-specialist can interpret the prediction map and make decisions about where an unknown material may originate from within Aotearoa using their own criteria. The main limitation with applying the isoscape to provenancing investigations is that the predictions highlight many potential regions-of-origin that have similar 87 Sr/ 86 Sr values when using the top 33% threshold. The threshold can be set to 20% or 10% but, as demonstrated here, these thresholds reduce the accuracy of origin predictions. Instead, predicted areas could be further constrained by combining the predictive strength of 87 Sr/ 86 Sr with other isotope systems, such as δ 2 H and δ 18 O [37, 74, 89] . Isoscapes for δ 2 H and δ 18 O already exist for Aotearoa and can easily be combined in a dual-isotope approach to predict the region-of-origin for unknown materials using the built-in assignment features of the assignR package [104, [136] [137] [138] . As people, animals, and materials are transported across increasingly large distances in a globalized world, issues with biosecurity and food security are rising. The global pandemic of COVID-19 has highlighted the difficulties we face in tracing and controlling the origins of animals and products that can transport viral loads. In Aotearoa, there are many endemic species and a strong local agricultural industry that must be protected from biosecurity threats. There is an urgent need to have tools which enable the provenancing of pests and agricultural products arriving in and being transported around the islands of Aotearoa, as well as confirm origins of products advertised as 'New Zealand made' that enter into overseas markets to protect valuable commodities from fraud. Isotopes are ubiquitous markers of provenance that are increasingly used to trace the origin of food or animals. In this study, we introduced the first bioavailable 87 Sr/ 86 Sr isoscape for Aotearoa and demonstrated how the isoscape can be used to certify the origin of agricultural products. We improved upon existing methodology to develop a bioavailable 87 Sr/ 86 Sr isoscape using the best available geospatial datasets to tune the regional isoscape. As anticipated, the primary drivers of bioavailable 87 Sr/ 86 Sr variability are the underlying geology, soil pH, and aeolian (dust and sea salt) deposits. We then tested and proved that there is potential for utilizing 87 Sr/ 86 Sr isotopes to determine the origin of cow milk and other agricultural products in Aotearoa. Currently, there is little to no geo-referenced bioavailable 87 Sr/ 86 Sr data derived from plants, animals, or soil in Aotearoa, beyond this study. As more data become available, the bioavailable 87 Sr/ 86 Sr isoscape model can be recalibrated and further improved, though we do acknowledge that sampling to fill current geographic gaps in the model is a costly undertaking. With the availability of this baseline bioavailable 87 Sr/ 86 Sr isoscape, our hope is to promote further research using 87 Sr/ 86 Sr isotopes in Aotearoa and abroad. Current provenancing methods in Aotearoa rely heavily on trace elements and light stable isotopes (δ 2 H, δ 15 N, δ 18 O, δ 13 C) to predict region-of-origin. The information provided by light isotopes and trace elements are useful for excluding potential regions of origin from consideration, but predictions can be further constrained if combined with other geologically derived isotope systems, like 87 Sr/ 86 Sr [37, 89] . Future provenancing projects in Aotearoa should implement a multi-isotope approach that uses 87 Sr/ 86 Sr alongside other light stable isotopes to produce both accurate and precise region-of-origin predictions. Supporting information S1 File. Sample preparation and analysis. S1.1. Aotearoa Plant Sample Preparation. S1.2. Aotearoa Topsoil Sample Preparation. S1.3. Collected Plant and Topsoil MC-ICP-MS Analysis. S1.4. Additional Plant and Topsoil Preparation and Analysis. S1. 5 Isotopes and trace elements as natal origin markers of Helicoverpa armigera-an experimental model for biosecurity pests Natal origin of the invasive biosecurity pest, brown marmorated stink bug (Halyomorpha halys: Penatomidae), determined by dual-element stable isotope-ratio mass spectrometry Plant invasions in New Zealand: global lessons in prevention, eradication and control Predator-free New Zealand: conservation country Human dimensions in the management of invasive species in New Zealand Can spectroscopy geographically classify Sauvignon Blanc wines from Australia and New Zealand? Biosecurity Pests and diseases we want to keep out of New Zealand. Ministry of Primary Industries Counterfeiting in the Primary Industry Sector and the threat to New Zealand's economy The Estimates of Appropriations for the Government of New Zealand-Primary Sector. Treasury. 2020 Correlation of geographical location with stable isotope values of hydrogen and carbon of fatty acids from New Zealand milk and bulk milk powder Identification of goat milk powder by manufacturer using multiple chemical parameters Regional classification of New Zealand red wines using inductively-coupled plasma-mass spectrometry (ICP-MS) Testing for geographic origin of fetal bovine serum Food traceability using the 87 Sr/ 86 Sr isotopic ratio mass spectrometry Applying the principles of isotope analysis in plant and animal ecology to forensic science in the Americas Bone deep: variation in stable isotope ratios and histomorphometric measurements of bone remodelling within adult humans Collagen turnover in the adult femoral mid-shaft: Modeled from anthropogenic radiocarbon tracer measurements Collagen turnover and isotopic records in cortical bone Analysing Sr isotopes in low-Sr samples such as single insects with inductively coupled plasma tandem mass spectrometry using N 2 O as a reaction gas for in-line Rb separation Establishing a strontium isotope baseline in New Zealand for future archaeological migration studies: A case study Strontium isotope composition of skeletal material can determine the birthplace and geographic mobility of humans and animals Statistical and geostatistical mapping of precipitation water isotope ratios Geographic assignment with stable isotopes in IsoMAP A simplified GIS approach to modeling global leaf water isoscapes Isoscapes: Understanding Movement, Pattern, and Process on Earth through Isotope Mapping Approaches to plant hydrogen and oxygen isoscapes generation Using isoscapes to model probability surfaces for determining geographic origins Determining geographic patterns of migration and dispersal using stable isotopes in keratins Mapping 87 Sr/ 86 Sr variations in bedrock and water for large scale provenance studies Application of stable isotopes and geostatistics to infer region of geographic origin for deceased unidentified Latin American migrants Spatial variation in bioavailable strontium isotope ratios ( 87 Sr/ 86 Sr) in Kenya and northern Tanzania: Implications for ecology, paleoanthropology, and archaeology Sr/ 86 Sr in European soils: A baseline for provenancing studies A strontium isoscape of north-east Australia for human provenance and repatriation Mapping multiple source effects on the strontium isotopic signatures of ecosystems from the circum-Caribbean region The first large-scale bioavailable Sr isotope map of China and its implication for provenance studies Strontium isoscapes in The Netherlands. Spatial variations in 87 Sr/ 86 Sr as a proxy for palaeomobility Mapping of bioavailable strontium isotope ratios in France for archaeological provenance studies Isotope values of the bioavailable strontium in inland southwestern Sweden-A baseline for mobility studies Sr/ 86 Sr maps as targeted strategy to support wine quality Determination of the source of bioavailable Sr using 87 Sr/ 86 Sr tracers: a case study of hot pepper and rice Strontium isotope analysis on cremated human remains from Stonehenge support links with west Wales Towards a biologically available strontium isotope baseline for Ireland A geostatistical framework for predicting variations in strontium concentrations and isotope ratios in Alaskan rivers A bioavailable strontium isoscape for Western Europe: A machine learning approach Triple sulfur-oxygen-strontium isotopes probabilistic geographic assignment of archaeological remains using a novel sulfur isoscape of western Europe Mapping and defining sources of variability in bioavailable strontium isotope ratios in the Eastern Mediterranean Variations of bioavailable Sr concentration and 87 Sr/ 86 Sr ratio in boreal forest ecosystems: Role of biocycling, mineral weathering and depth of root uptake 87 Sr/ 86 Sr ratios in modern and fossil food-webs of the Sterkfontein Valley: implications for early hominid habitat preference Geochemical baseline soil surveys for understanding element and isotope variation across New Zealand Soil Quality-Extraction of Trace Elements From Soil Using Ammonium Nitrate Solution Predictors of pesticide concentrations in freshwater trout-The role of life history Prehistoric migration at Nebira, South Coast of Papua New Guinea: New insights into interaction using isotope and trace element concentration analyses Origin Authentication of Bovine Milk Constraints on assembly of Tongariro and Ruapehu andesite magmas based on Sr-isotope compositions of plagioclase and groundmass Does the trace element composition of brown trout Salmo trutta eggs remain unchanged in spawning redds? Evaluation of a strontium-specific extraction chromatographic method for isotopic analysis in geological materials Isotopes and trace elements as geographic origin markers for biosecurity pests Decision intelligence: determining pest natal origin. Final report for contract PBCRC2111 Plant Biosecurity Cooperative Research Centre 37 Infer Geographic Origin from Isotopic Data. R package version 1.1.3 A test of geographic assignment using isotope tracers in feathers of known origin Contrasting assignment of migratory organisms to geographic origins using long-term versus year-specific precipitation isotope maps Random forests New Zealand Geology: an illustrated guide Strontium isotope geology Late Quaternary slip rates and slip partitioning on the Alpine Fault, New Zealand The geology of the West Coast from Abut Head to Milford Sound-Part 1 Overview of the Median Batholith, New Zealand: a new interpretation of the geology of the Median Tectonic Zone and adjacent rocks New Zealand's geological foundations Geochronology and geochemistry of a Mesozoic magmatic arc system The isotope and trace element budget of the Cambrian Devil River Arc System, New Zealand: identification of four source components U-Pb geochronology of mid-Paleozoic plutonism in western New Zealand: Implications for S-type granite generation and growth of the east Gondwana margin Geological map of New Zealand 1:1000000. GNS Science Geological Map 2. 2 sheets. Lower Hutt: GNS Science New Zealand: redefined Major periods of erosion and alluvial sedimentation in New Zealand during the Late Holocene Stratigraphy, age, and correlation of voluminous debris-avalanche events from an ancestral Egmont Volcano: implications for coastal plain construction and regional hazard assessment Geology of the Aoraki area: scale 1:250,000. Institute of Geological & Nuclear Sciences 1:250,000 geological map 15, 71 p. + 1 folded map. Lower Hutt: GNS Science Structure of the alpine schists of South Westland Ar-Ar dating of K-feldspar in low grade metamorphic rocks: Example of an exhumed Mesozoic accretionary wedge and forearc Geochemical characterization, evolution and source of a Mesozoic accretionary wedge: the Torlesse terrane Overview of New Zealand's Climate The climate and weather of West Coast Quantifying the economic cost of invertebrate pests to New Zealand's pastoral industry Classifying the introduction pathways of alien species: are we moving in the right direction? Corruption, development and governance indicators predict invasive species risk from trade Import volumes and biosecurity interventions shape the arrival rate of fungal pathogens The accidental introduction of invasive animals as hitchhikers through inanimate pathways: a New Zealand perspective The stable isotope ratios of marijuana. II. Strontium isotopes relate to geographic origin Profiling of heroin and assignment of provenance by 87 Sr/ 86 Sr isotope ratio analysis Multi-isotope approaches for region-of-origin predictions of undocumented border crossers from the US-Mexico border: Biocultural perspectives on diet and travel history Investigating human geographic origins using dual-isotope ( 87 Sr/ 86 Sr, δ 18 O) assignment approaches Precipitation isoscapes for New Zealand: enhanced temporal detail using precipitation-weighted daily climatology A δ 2 H Isoscape of blackberry as an example application for determining the geographic origins of plant materials in New Zealand A feather-precipitation hydrogen isoscape model for New Zealand: implications for eco-forensics A special thanks to the Centre for Trace Elements team at the University of Otago for the countless hours of lab guidance and supervision. Additional thanks to AsureQuality for assistance collecting B3 samples and to David Murphy, Queensland University of Technology for managing analysis of the B3 samples. We acknowledge the SPATIAL (Spatio-temporal Isotope Analytics Lab, University of Utah) Short Course for providing the environment to design this work and the funds to carry it out to RTK.