A strontium isoscape of Italy for provenance studies

1. Department of Cultural Heritage, University of Bologna, Via degli Ariani 1, 48121 Ravenna, Italy. 3 2. Department of Chemical and Geological Sciences, University of Modena and Reggio Emilia, 41125 Modena, Italy. 4 3. Lamont-Doherty Earth Observatory, Columbia University, Palisades, New York 10964, USA. 5 4. Dipartimento di Storia Culture Civiltà, University of Bologna, 40124 Bologna, Italia. 6 5. Department of Archaeology, Durham University, Durham, United Kingdom. 7 6. Max Planck Institute for Evolutionary Anthropology, Department of Human Evolution, 04103 Leipzig, Germany. 8


29
Geology is biological destiny: Whatever minerals land or are deposited in a place determine what or who can make 30 a living there millions of years later.
('bioavailable' + 'non-bioavailable'; see Table S1). This allowed us to obtain a broader overview of the Sr isotope 110 distribution across Italy. 111 Novel data (n = 89) were generated from modern environmental and archaeological samples by solution MC-113 ICPMS analyses. Samples include modern vegetation, archaeological and modern teeth, snails, waters, rocks and 114 soils. These samples are from areas where archaeological studies are in progress and thus were integrated into the 115 database. Five meteoric water samples collected from pluviometers located in the Emilian Apennine 116 (Montecagno, 44°19'57.76'' N; 10°21'58.57'' E) were also measured for their Sr isotopic composition. These values 117 were not included in the spatial model, but are presented as possible end-members for the Sr cycle in the biosphere, 118 possibly helpful for future studies on Sr mixing (Table 2). 119

Solution MC-ICPMS 121
Samples were processed at the Geochemistry Lab of the Department of Chemical and Geological Sciences 122 (University of Modena and Reggio Emilia) and at Durham University. All the reagents employed were of suprapur 123 grade (Modena) or bidistilled (Durham). Biominerals (i.e. teeth and snail shells) were cleaned with MilliQ water 124 and digested using concentrated HNO3. The bioavailable Sr fraction from soils instead was extracted using 0.25M 125 acetic acid. Bulk rocks samples were totally digested using a mixture of concentrated HNO3 and HF. Waters were 126 filtered and acidified with HNO3 to a concentration of 3M. After drying and re-dissolution by 3M HNO3, all 127 samples were processed using the Eichrom Sr-spec resin. The 87 Sr/ 86 Sr ratios were determined by Neptune  Emilia and one at the Northern Centre for Isotopic and Elemental Tracing at Durham University. Detailed  where the isotope groups were defined using clustering techniques on the data itself, we relied on a simplified 157 geological map of Italy (Figure 1 high number of isotope data is available in the literature for metamorphic and magmatic rocks across Italy, which 165 have been measured to understand the geodynamic events that led to the formation of the Alps and Apennines 166 and their emplacement at crustal level. Although most of these data were not included in the database, because no 167 geolocalization was available, their isotope signature was used to define isoclasses as building blocks of the Italian 168 Sr isomap. In addition, several published Sr isotope data were measured on single mineral phases and therefore, 169 being not representative of the bulk rock, could not be used for our purpose.  formations. Isoclass 7 (expected 87 Sr/ 86 Sr > 0.70920) includes old metamorphic and magmatic rocks of the 178 crystalline basement and younger volcanics whose magmatism is affected by a radiogenic Sr isotope source. Isoclass 179 8 finally includes all the geolithologies that we were not able to attribute to an isotope class due to their hybrid 180 nature (i.e. siliciclastic rocks) or to their wide Sr isotope variability (i.e., Permian to Devonian carbonates have a 181 very wide range of Sr isotope ratios across several of our defined classes). 182 In attributing the isoclass to a particular geolithology or formation we confronted local rock values from literature 183 and, whenever possible, double checked their consistency with the bioavailable values of our database. When no 184 data were available, we considered the type of rock (i.e. mineralogy) and the age of formation. Initially, we defined 185 several more isoclasses in the Sr isotope range especially in the range between 0.7092 and very radiogenic values 186 (up to 0.75). However, we could attribute with certainty only a few data points from Sardinia to these classes, and 187 therefore we finally grouped all Sr isotope ratios > 0.7092 in a unique class (isoclass 7). We stress that the 188 attribution of an isoclass has not been arbitrary and any attribution is either backed up by isotopic data or 189 consistent with a particular type of magmatism or deposition event (e.i. seawater curve for marine carbonates of cyclical-like structure, with a first maximum located at approximately 250 km ( Figure S2). The prediction power 194 of the models was evaluated using a 10-fold cross-validation method through SAGA 7.9 (Table 3). The 195 interpolated Kriging models were imported into QGIS 3.18 to generate the final distribution maps (freely available 196 online at geochem.unimore.it/sr-isoscape-of-italy). We note here that Sardinia was excluded from the Ordinary 197 Kriging due to the low number of data from the area. 198

Data description and distribution 201
Descriptive statistics for the data considered in this study are reported in Tables 1 and S1 and summarized in 202  Table 1). Bioavailable samples show an average 87 Sr/ 86 Sr ratio of 0.70941 ± 0.00632 (2 SD), 207 and span between 0.70354 and 0.76384, with a median value of 0.70883. The kernel density distribution of the 208 bioavailable data is strongly asymmetric and leptokurtic (skewness = 8.14; kurtosis = 99.16). Notably, the non-209 bioavailable samples, including all the rocks and bulk soils, display an average 87 Sr/ 86 Sr ratio of 0.71069 ± 0.01054 210 (2 SD), ranging between 0.70319 and 0.75300, with a median value of 0.70900 (Table S1). The distribution of 211 the non-bioavailable dataset is asymmetric but less leptokurtic than the bioavailable (skewness = 3.93; kurtosis = 212 22.40). Yet, we stress that the number of non-bioavailable data (n = 352) here considered is remarkably lower than 213 the data in the bioavailable dataset (n = 1568), potentially influencing our observations on the data. Similarly, the 214 uneven spatial distribution of 'non-bioavailable' samples across Italy certainly influenced data evaluations and use 215 for this class. 87 Sr/ 86 Sr ratios of the bioavailable samples were also exploratively plotted against latitude and 216 longitude ( Figure 3), searching for potential correlations between these variables. However, no statistically 217 significant trend was observed (both R 2 < 0.1). Yet, the two graphs clearly show a preferential distribution of the 218 highest radiogenic Sr values northwards (latitude 44-47° N) and eastwards (longitude 7-12° E). This is expected 219 due to the presence of old metamorphic and magmatic rocks in the Alpine area and magmatic-metamorphic 220 provinces in Central Italy (Tuscany, Latium), and also evident when data are plotted by Italian macroregions 221 ( Figure 2). 222 Five meteoric waters, not included in the previous statistics evaluations (and the interpolated maps) range between 223 0.70848 and 0.70924, and represent an end-member of the Sr bioavailable cycle. These five waters were sampled 224 from the same pluviometer located in the Emilian Apennine, and they were seasonally collected ca. 3-to-5 months 225 apart from each other. These data highlight a remarkable temporal variability of the local rainwater likely due to the changing contribution of seawater aerosol and crustal dust, with a possible important influence on the local 227 bioavailable Sr (Négrel et al., 2007).

Maps 242
The isoclass map of Italy ( Figure 1) allows a first order distinction between the radiogenic Sr isotope provinces, 243 related to the 'old' crustal and radiogenic Sr isotope magmatism units mainly present in the Alps, Calabria, 244 Sardinia and Central Italy, and the unradiogenic provinces related to the depleted mantle magmatism mainly in 245 the Southern Alps and Sicily. Yet, more information can be gathered through the isoscape maps (Figures 4 and 5). 246 These were built modelling the two datasets, namely 'bioavailable' and 'all'. Each figure includes two maps 247 obtained with two distinct Kriging approaches: Ordinary and Universal with external drift. The evaluation of 248 performance of the two models is reported in Table 3. Both methods produced satisfying results, with relatively 249 low normalized root mean squared errors (NRMSE ~3-4%), explaining between ~60 and ~70% of the isoscape 250 variance (R 2 ). In general, Universal Kriging (with eternal drift) seems to outperform Ordinary Kriging, although 251 the difference is not remarkable ( Table 3). The lowest RMSE is observed for the 'bioavailable' Universal Kriging, 252 and is equal to 0.0020; instead, the highest RMSE (0.0024) was obtained for the 'all' Ordinary Kriging model. prediction power of the Kriging method, both in terms of data over-fitting (higher R 2 ) and worse variogram 255 modelling (see also Figure S4). To further evaluate the prediction of our modelling we measured the prediction 256 standard errors for the Kriging maps ( Figure S4). Both models (i.e. Ordinary and Universal) show similar standard 257 prediction errors, ranging from ca. 5E-7 to 5E-6 for the 'bioavailable' dataset and from 2E-7 to 2E-5 for the 'all'   The largest differences in terms of isoscape predicted values among the 'all' and the 'bioavailable' maps arise indeed 277 in these areas (particularly Tuscany and Latium), due to the presence of even higher radiogenic values in local 278 rocks, only partially identified in the bioavailable pool (see Figure S5). The north-western Alpine area also shows 279 significant differences (both in negative and positive) between the two datasets. However, here, only few rock western Italy) are probably linked to model's predictions inaccuracies rather than actual variations of the 87 Sr/ 86 Sr 282

ratio. 283
Overall, several small 'hotspots' (both negative and positive) can be recognized when comparing the predictions 284 of the two datasets, particularly in the Alps. We stress that the number of samples in these areas is lower than in 285 other localities; however, another explanation might lie in the complex geometry of the Alps where the 286 bioavailable Sr isotope ratios might differ from those of the exposed rocks because of the geological complexity of 287 the nappes that overthrust each other in the belt and therefore in the differential contributions to the bioavailable 288 Sr possibly from other reservoirs. 289 Sharper details of the isotope zones can be observed in the Universal Kriging map compared to the Ordinary 290 Kriging, due to the definite isoclass boundaries of the guiding map. In general, when looking at specific areas of 291 the map, the Universal Kriging model should be more accurate in terms of spatial prediction, particularly for those 292 areas with few data available. However, the Ordinary Kriging map seems to better mimic the natural averaging of 293 Sr isotope values due to weathering and mixing processes.

346
Overall, these data suggest that soils (leachates) and plants best reflect the local bioavailable Sr pool, although 347 possibly contaminated by modern and/or antropic end-members. Fauna enamel, if truly local as in the case of 348 domesticated macro-mammals or small home range micro-mammals, mixes various bioavailable Sr sources and 349 more closely mimics the local food and drinking sources. Such evidence clearly highlights the intrinsic limits in 350 using isoscapes, which are commonly composed by a patchwork of literature data from different samples, or 351 modelled on specific samples collected ad hoc (as soils or plants). Yet, we stress here that Sr isotopes need to be et al., 2021). This, in turn, suggests that provenancing through isoscapes, and isotope baselines in general, need to 354 be performed with caution. Hence, isoscapes must be considered as 'guides' for data interpretation, rather than 355 an unequivocal provenancing tool, justifying their composite nature to better understand the variability of local 356 We took advantage of the produced maps to discuss a local case study and the definition of local baseline in 367 archaeological studies, a currently hot-topic within the field of provenance and mobility studies. In this sense, 368 regional and (extra)national isoscapes are key in understanding the variability of the local Sr pool, broadening our 369 understanding on the mixing of the different end-members to obtain certain isotope signatures in (geo)biological 370

samples. 371
Distribution maps of Sr isotopes provide a solid interpretative basis for provenance and traceability studies. Our 372 maps and database are freely accessible online and will be updated in the future when new data become available. 373 In this sense, we will continue to collect and analyse new environmental samples from low-density areas (such as 374 Sicily and Sardinia) to improve the prediction power of the models. In addition, we plan to employ novel methods 375 for the spatial modelling of isotope data, using different predictors and machine learning approaches.   Figure S3. Box plot graph representing the 87 Sr/ 86 Sr ratios of the different sample categories. The gray shadows are data distributions obtained by kernel density estimations. Note that the graph is cut at 0.724 to improve readability, with some outliers plotting beyond the graph range. Figure S4. Prediction standard error for the Kriging models. Maps were obtained using SAGA 7.9 and QGIS 3.8. Figure S5. 87 Sr/ 86 Sr difference between 'all' and 'bioavailable' Universal Kriging models; e.g. a positive value means that Sr isotope data are higher in the 'all' map. Maps were obtained using SAGA 7.9 and QGIS 3.8.