Tooth enamel oxygen “isoscapes” show a high degree of human mobility in prehistoric Britain

A geostatistical model to predict human skeletal oxygen isotope values (δ18Op) in Britain is presented here based on a new dataset of Chalcolithic and Early Bronze Age human teeth. The spatial statistics which underpin this model allow the identification of individuals interpreted as ‘non-local’ to the areas where they were buried (spatial outliers). A marked variation in δ18Op is observed in several areas, including the Stonehenge region, the Peak District, and the Yorkshire Wolds, suggesting a high degree of human mobility. These areas, rich in funerary and ceremonial monuments, may have formed focal points for people, some of whom would have travelled long distances, ultimately being buried there. The dataset and model represent a baseline for future archaeological studies, avoiding the complex conversions from skeletal to water δ18O values–a process known to be problematic.

Scientific RepoRts | 6:34986 | DOI: 10.1038/srep34986 follows that the bioapatite isotopic composition is ultimately correlated with that of local environmental waters. Equally, soft tissue δ 18 O values are correlated with those of local waters.
Ehleringher et al. 2 demonstrate that mobility of modern individuals can be traced by means of the δ 18 O and δ 2 H of hair keratin. These authors provide a useful theoretical and visual method to predict δ 2 H and δ 18 O values of human hair by mapping the expected values in a Geographical Information System (GIS) platform. Their predictive model is based on expected tap water values from across different areas of the USA and the mass balance between the dietary oxygen and hydrogen isotope inputs from a general "supermarket diet", i.e. a modern globalised diet. In a comparable approach, we have used the isotopic characteristics of skeletal remains to identify human mobility of past human populations by mapping the variation in tooth enamel δ 18 O p of archaeological remains coming from across Britain. The resultant map is based on the expectation that prehistoric people would have sourced their dietary and drink supplies locally, and that these sources are compositionally related to the mean isotopic values of rain and groundwater. This "isoscapes" (isotopes in the landscape) approach is novel in archaeology.

Environmental-based Approaches to Mobility
Many isotopic studies that have employed δ 18 O p to trace human mobility in archaeological or forensic contexts have focused on measuring the values of single individuals to compare them with the isotope ranges of local waters expected at the site of recovery 12,[17][18][19][20][21] . When such isotopic values differ from the ones expected, individuals have been considered 'non-locals' to the area 13 . Some studies have further attempted to source individuals to specific places of origin, or places where they would have spent their childhood (i.e. the period during their tooth mineralization). However, one of the main difficulties with this process is discerning between the many areas with the same combination of isotopes across the landscape, both within any one country, and across different countries 22 .
Standard practice for investigating mobility of past people with isotopes is to measure the isotopic composition in the body tissues (for example 18 O/ 16 O) and to compare this composition with the distribution of the same isotopes in the local and wider environment. In most archaeological studies the material chosen for analysis is the skeleton, which is the tissue most prone to survive in the burial environment. Bioapatite phosphate and carbonate oxygen isotopes (δ 18 O p or δ 18 O c ), measured in teeth or bones, are used to trace human mobility, by converting these values to their equivalent water δ 18 O w values according to specific conversion equations available in the literature 15,16,[23][24][25] . These equations take into account the physiological effects during the integration of the different isotopes into the tissues. The "converted" values so calculated are subsequently compared to rainwater or groundwater patterns of variation, either obtained from the literature [26][27][28][29] or extrapolated from the general distributions of isotope ratios in precipitations 30 . With this method, non-local individuals are identified based on the difference of their skeletal δ 18 O values (converted to water values) from that of the drinking water in their place of recovery. It must be considered, however, that the isotope composition of the same elements in the environment and in living individuals does not vary in exactly the same way. The relationship between the measured value in an individual and the value expected on the basis of the place of recovery is therefore not straightforward, and the biological values recorded in local populations may be more or less consistent with those predicted. Apart from the possibility that someone in the group is not local to the area, the variability in oxygen isotopes recorded in people living in a specific area could be due to several other factors: (1) Some individuals may have been affected by short-term climate conditions (warmer/colder, wetter/drier periods) occurring during childhood formation of their teeth. This may lead to atypical δ 18 O values ( 18 O-enriched or 18 O-depleted). Mean annual water values, with which these are compared, are averaged over a long period of time (normally 10 to 30 years), a period which is longer than that required for the tooth to mineralise; (2) Sourcing drinking water from reservoirs other than the local groundwater, for example from rivers coming from higher latitudes or from lakes or ponds, may also contribute to altering the individuals' expected skeletal δ 18  These factors can all contribute to altering the direct relationship between individuals' oxygen isotope ratios and the environmental ratios of their place of origin. The best approach is therefore to avoid the conversion of skeletal δ 18 O p to water δ 18 O w , and instead compare the skeletal values directly with other human phosphate values.
This paper provides a distribution model for human tooth enamel δ 18 O p values across Britain and describes its most likely variation during the Chalcolithic-Early Bronze Age. The map of variation produced is a predictive model of skeletal δ 18 O p values across the region and will allow a direct evaluation of δ 18 O p measured in individuals of comparable age, with no need to convert their δ 18 O p values to water δ 18 O w values. The model assumes that the majority of the samples investigated represent individuals that were "local" to the place of recovery. It also assumes that sampled humans did not engage much in dietary behaviours such as drinking stewed or brewed beverages, which would modify substantially their δ 18 O p values or, if they did, they would do similarly at the population level. Finally, it is assumed that environmental conditions remained consistent over the period of the Chalcolithic and Early Bronze Age.

A Geostatistical Approach to Mobility
The identification of human mobility based on isotope compositions is contingent upon characterisation of the spatial variation of the expected isotope values. Without adequate characterisation of this spatial dimension, it is not possible to identify individuals who are 'local' to an area and those who are not with any degree of certainty. Teeth from 261 individuals (denoted by the prefix SK) were sampled in order to characterise spatial variation in human skeletal δ 18 O p for the Chalcolithic and Early Bronze Age (2500-1500 BC) as part of the Beaker People Project 7 . The individuals analysed include 22 samples that fall outside Chalcolithic-Early Bronze Age period (see Supplementary Table S1). All are from Britain except for SK 315 (Burial 30, Mound of Hostages, Tara). This, the only sample from the Republic of Ireland, has been described in detail elsewhere 32 and is excluded from the modelling below. The sample locations indicated in Fig. 1 reflect a number of geographic biases that must be acknowledged prior to undertaking any form of spatial analysis: (1) Sample locations correspond to the places where individuals were buried rather than the places where they lived: in some cases, Chalcolithic and Early Bronze Age settlements may have been located several kilometres from the place of burial 33 ; (2) Many of the skeletal samples were recovered from burial mounds excavated by nineteenth-century antiquaries, the locations of which were sometimes poorly recorded. Sample locations were recorded as 4-figure National Grid References, i.e. they were located to the nearest 1km grid square-the highest level of precision common to all samples; (3) Most samples were recovered from areas with geology favourable to the preservation of human bones (e.g. chalk, limestone). Skeletal remains tend not to survive in acidic soils such as those in much of Wales and Scotland. The majority of the Scottish samples included in this data set survived due to the burial rite, which involved a stone cist and avoided soil contact 7 . (4) The resultant dataset exhibits preferential sampling due to the survival of skeletal remains, the availability of material from museum collections, and the areas in which excavators and antiquaries were active. A large number of samples are from the south of England and the northeast of Scotland, with very few samples from the north and centre of Wales and the south and west of Scotland.
Given the complexity of the relationship between δ 18 O p and δ 18 O w , the identification of individuals who are 'local' or 'non-local' to a particular area is dependent here upon the comparison of the δ 18 O p values with all the  18 O p values for each of these locations was more than half of that for the whole of Britain. Despite this variability-a defining characteristic of the spatial variation in δ 18 O p values that should not be ignored or played down-the data exhibit positive spatial autocorrelation at the 99% confidence level. The Moran's I test statistic, a global measure of spatial autocorrelation, was calculated for the δ 18 O p values for the individuals sampled as part of the Beaker People Project (I = 0.315) and compared to the expected value (I E = − 0.004). The z-score (z = 6.004) and p-value (p = 0.000) indicate that the test statistic is significant at the 99% confidence level. Spatial autocorrelation, the assumption that the values of nearby samples are closely related to one another 34 , is a pre-requisite of the geostatistical techniques that underpin the approach used in this paper.
Cluster/outlier analysis was carried out in order to characterise spatial patterning in the δ 18 O p values. Anselin's Local Moran's I test statistics were calculated for each of the samples, with positive values of the test statistic indicating similarity and negative values indicating dissimilarity 35 . Spatial relationships between neighbouring features were defined by inverse distance, with a threshold of c.85km (the distance at which each sample had at least one neighbour). Row standardization was applied to the spatial weights used in the calculation of the test statistic for each sample. Corrections based on the False Discovery Rate were applied to the classification of samples corresponding to clusters/outliers 36 . The values of the test statistic are used to identify samples that form part of a cluster of high/low values (HH/LL) and samples that correspond to outliers with unusually high/low values (HL/LH). The cluster/outlier analysis (see Supplementary Fig. S1) reveals an underlying trend in the δ 18 O p values, with clusters of high values in the south and west of Britain (Wales and South West) and clusters of low values in the northeast of Britain (Eastern Scotland and North Eastern Scotland). Two outliers, samples with unusually high or low δ 18 O p values, were also identified; namely the samples from Woodhenge (South West) SK308 and 309. These samples, which could be interpreted as samples from 'non-local' individuals, were excluded from further analysis.
A continuous surface predicting the spatial variation in δ 18 O p values (Fig. 1a) was subsequently generated using Ordinary Kriging-a geostatistical interpolation technique that estimates values for non-sampled locations based on weighted averages of the values of nearby samples [37][38][39] . Sample locations were de-clustered in order to compensate for any biases introduced by preferential sampling. The underlying trend identified through cluster/outlier analysis was modelled using a 2 nd order polynomial trend surface. This trend was subtracted from the data, with the semi-variance modelled on the residuals, and added back in at the end of the interpolation process. The semi-variance was modelled using a stable model (nugget = 0.299684, partial sill = 0.015628, range = 121,384.818021 and parameter = 1.170313), with a smoothed search neighbourhood (smoothing  26 ).
Scientific RepoRts | 6:34986 | DOI: 10.1038/srep34986 factor = 0.5). Cross-validation of the model, with samples excluded on a case-by-case basis and the observed values compared to the value predicted from the remaining samples, indicates that the model is accurate and unbiased. The standard error surface derived from cross-validation (Fig. 1b) indicates that the predicted values are typically a good fit to the measured values in locations from which multiple samples were taken.
Cluster/Outlier analysis and Ordinary Kriging were carried out using ArcGIS 10.3.1. The parameters used to model the trend and semi-variance were optimised using the Geostatistical Wizard which forms part of the Geostatistical Analyst extension. Geostatistical Analyst supports Ordinary Kriging using either a constant mean or a deterministic trend.

Spatial Variation in δ 18 O p values
The δ 18 ± 0.7). The two datasets are statistically indistinguishable. The minor difference at the higher end of the ranges could be due to the larger sample size in this study or differences in the spatial distributions of samples from the respective studies.
The underlying trend in the point data, with high values in the south and west of Britain and low values in the northeast of Britain, is clearly evident within the prediction surface (Fig. 1). This trend is broadly in line with that recorded for the rain and groundwater oxygen isotope ratio distribution for Britain and Ireland 26,41 . In Britain, rain and groundwater δ 18 O values decrease from south to north and from west to east, due to the continuous depletions of air masses as they move over the region. A digitised version of the groundwater isotope map is provided in Fig. 2. The gradient of depletion in 18 O is quite steep in the western regions, where there is also the highest amount of rainfall, and more gentle towards the eastern regions. Provided that the climate conditions and air mass circulations have not changed substantially in the Holocene, as suggested by Darling et al. 26 , this map should provide a good representation of the natural waters available to individuals who inhabited these areas. As stated, individuals who lived in one place for their entire life and who obtained water and food from local sources should have tissue isotope ratios similar to that of the local environmental water (plus or minus a physiological fractionation factor, which is accounted for in the conversion equations 15,16,[23][24][25]. Theoretically, if all of the individuals sampled for the Beaker People Project had been local, and unaffected by any of the factors described above, the geographic trend in the distribution of skeletal δ 18 O p values should mirror that of the water (δ 18 O w ). Although the general trend of depletion in δ 18 O values observed for the natural waters (Fig. 2) broadly persists in the skeletal isotopes (Fig. 1a), key differences can be seen between the modelled spatial variation of δ 18  These key differences can in part be attributed to factors that may affect the relationship between isotope ratios in living populations and the environment, as well as biases in the spatial distribution of samples discussed above. However, the overall variation in the spatial distribution of δ 18 O p suggests that something else besides the physiological/environmental offsets between the local water and the skeletal isotopic values is playing a role.
Already, the two individuals with unusually high (δ 18  Although it is tempting to infer that these individuals were immigrants who brought Beaker pottery to Britain from continental Europe, any notion of a mass migration of the Beaker people is overly simplistic. In fact, the Woodhenge burials date to long after the 'first-wave' arrival of Beakers; SK 272 from Staxton also dates to later than the start of the Beaker package 7 and other burials did not contain this form of pottery. Instead of representing 'first-wave' migrants from continental Europe, the outliers and isolated clusters can be considered to reflect a high degree of mobility during the Chalcolithic and Early Bronze Age-with movement both within Britain and between Britain and the Continent. This helps to place in context the apparent long-distance mobility of the so-called Amesbury Archer, a Beaker individual previously identified on the basis of his δ 18 O p value as a continental migrant to Britain, although not quite early enough to have arrived in a first wave 42 . The results of the Beaker People Project indicate that his long-distance mobility was exceptional rather than typical of these patterns of movement 6 . The underlying trend in predicted values and small standard prediction errors (Fig. 1) masks the variability identified in the measured values of samples from the same region. Plotting the residuals from cross-validation of the geostatistical model, i.e. the differences between the predicted and measured values for each sample (Fig. 3), highlights the magnitude of this inter-regional variability. Samples with large residuals ( 18 O-depleted < − 0.6‰ and > 18 O-enriched 0.6‰, equivalent to over 2σ of the measurement error) are found in many areas, including the Stonehenge region (South West), the Peak District (East Midlands and West Midlands) and the Yorkshire Wolds (Yorkshire and the Humber). The large residuals for the samples from the Stonehenge region, the Peak District and the Yorkshire Wolds are comparable to those in samples recovered from Medieval ports and cities, based on forthcoming geostatistical modelling of data published by Evans et al. 40 . Each of the highlighted areas is rich in Chalcolithic and Early Bronze Age funerary and ceremonial monuments that would have formed focal points or gathering places for people, some travelling long distances and ending up being buried there.
The spatial outliers can therefore be regarded as non-local to the regions that they were found and, by extension, the high degree of variability identified in the Yorkshire Wolds, the Peak District and the Stonehenge region can perhaps best be accounted for by a high degree of population mobility in those areas.

Conclusion
This paper presents a geostatistical model of spatial variation in the distribution of human phosphate oxygen isotope values (δ 18 O p ) for Chalcolithic and Early Bronze Age Britain. The model can be employed to study the characteristics of local human populations and the mobility of individuals, and provides a baseline against which future samples from archaeological burials can be compared, without the need to undertake complex conversions of the data from skeletal to water δ 18 O values, a process known to be problematic. This model helped to identify δ 18 O p outliers, i.e. individuals that are statistically too different from the rest of the population to be 'local' , and can be used to determine whether measured values from new samples are 'local' or 'non-local' to the area from which they were recovered. The hypothesis underpinning this model is that the majority of people living in a certain area are not immigrants. It should be noted that the dataset used to create our geostatistical model is not complete and does not include data from all geographic regions. However, the dataset can be supplemented with samples from other regions, where available, once the offsets resulting from the lack of internationally-recognised materials for phosphate oxygen isotope analysis and commonly-agreed protocols for precipitation or calibration techniques are better understood. Although specifically created for Britain, the model can be exported to other settings or isotopes. The approach adopted in this study can be applied to datasets from other geographical regions or chronological periods. Based on oxygen isotopes alone (or the isotopes of any other individual element), it is not Scientific RepoRts | 6:34986 | DOI: 10.1038/srep34986 possible to distinguish, definitely, immigrants from local people if they have grown up in environments with similar isotope compositions to the place of recovery. It is therefore advisable to combine oxygen isotope profiles with other isotopic profiles, and/or other archaeological information, to be able to spot immigrants that can resemble locals based on one of their isotopic ratios. Other isotope measurements (Sr and S) are available for the individuals included in the Beaker People Project and are published separately 6,7,43 .

Methods and Materials
The model of δ 18 O p variation was created from a new database of individual teeth coming from Chalcolithic and Early Bronze Age human burials in Britain and Ireland. Most of this collection dates from 2500 to 1500 BC 6,7,44 and many are individuals that had been buried with Beakers, a type of pottery that appeared in Britain and Ireland shortly after 2500 BC 45 . A total number of 261 individual teeth were analysed. Tooth enamel fragments of about 20 mg were sampled from each tooth at the Division of Archaeology of the University of Bradford, prepared and precipitated as silver phosphate (Ag 3 PO 4 ) at the Department of Human Evolution-Max Planck Institute for Evolutionary Anthropology in Leipzig (DE), and measured for their δ 18 O p value back at the University of Bradford. Second permanent molar teeth were preferentially chosen although, in a few cases, when this tooth was not available, premolar or other teeth were used (SK51, SK53, SK129, SK135). Permanent molar teeth (except first molars) and premolars mineralize over a comparable period of time, approximately from 2-3 years of age to about 8, thus sensibly longer than one year 46,47 . The teeth chosen target a very early stage of the individual's life, i.e. childhood. The time span covered by the tooth mineralization should provide a good indication of the mean annual environmental conditions of the area where the individual grew up and take into account inter-annual climate variability. First molars were avoided for their likelihood of being characterized by breastfeeding isotopic signatures 48 . The sampling was standardised by selecting enamel sections from the same region in the tooth, i.e. the central part of the tooth crown. The investigated fragment should therefore have mineralised between the age of two and eight, and represents a weighted mean of the food and drink ingested over multiple years. Additional sampling and methods details can be found in Montgomery et al. 43 .
Tooth enamel bioapatite was precipitated as silver phosphate (Ag 3 PO 4 ) following a customised version of the procedures described in Dettman et al. 49 . Its isotopic composition was measured via Thermal Conversion Mass Spectrometry following the principles outlined in Gehre and Strauch 50 . Enamel chips were crushed to a homogenous fine powder in an agate mortar and pestle, and about 8-10 mg of this powder was dissolved in 2ml HF for 24h. Following the results in Grimes and Pellegrini 51 , no pre-treatment protocol was performed to clean the samples of potential organic contaminants. The solution containing dissolved phosphate ions was separated from the precipitated CaF 2 by centrifugation, and the solid CaF 2 rinsed three times with MilliQ water, the rinsing water being also added to the phosphate solution. The liquid was subsequently mixed with 2M NH 4 OH and AgNO 3 to precipitate Ag 3 PO 4 , which crystals were separated by centrifugation from the solution and dried out overnight at 60 °C. An aliquot of 250-350 μ g of these homogenised crystals was loaded into silver capsules and combusted in a Thermo Finnigan TM High Temperature Conversion Elemental Analyser, coupled with a stream of helium to a Thermo Finnigan TM Delta Plus XL mass spectrometer set at a temperature of 1450 °C.
The results given in the SI for each sample are the average of two to five repeated measurements. The average reproducibility of the isotopic measurements for each sample was 0.13‰ ± 0.09 (1σ ). Raw data were corrected with repeatedly measured standard materials, regularly analysed within the run in the same conditions as the samples, to compensate for drifts in both the mass spectrometer and the peripheral. Shifts from the assigned values of the standards were applied to the value of the measured sample 52 . The correction applied to phosphate oxygen isotope measurements is not free of complications, given the lack of internationally recognised standards 53,54 and the poor uniformity of the precipitation methods and calibration approaches employed in the different laboratories. For this project, raw data were corrected against widely-used silver phosphate materials (TU1 = 21.1‰ and TU2 = 5.5‰) that 'bracketed' the composition of all the samples investigated in this study (i.e. they have lower and higher values). These working standards were reported in Vennemann et al. 55 . A two-point linear correction based on these standards was applied, as also recommended by Meier-Augenstein 53 . A third reference material, Ag 3 PO 4 precipitated in large amounts in our lab from the phosphorite rock international standard SRM 120c (formally NBS 120c), was instead used as an external check for the calibration. Although frequently used as a reference standard for phosphate oxygen isotopes because it is a suitable inorganic analogue to bioapatite, this material is technically only a reference standard for its chemical content and not for its isotopic composition. Additionally, it has to be precipitated (i.e. manipulated) as Ag 3 PO 4 in each laboratory before isotopic analysis. It also has a quite large interval of reported values ranging from 19.9‰ to 22.6‰, with a mean of 21.7‰ ± 0.5‰ 20 . Finally, its rather 18 O-enriched value, outside the range of most phosphate values for humans from northern Europe and Britain, made it unsuitable to calibrate our dataset. Within this study, our SRM 120c Ag 3 PO 4 working standard gave a range of 22.1 ± 0.3‰ (n = 72), in line with reported values.
SRM 120c was also employed in this study to check the accuracy and reproducibility of the precipitation protocol by including samples of this standard in each precipitation batch together with the unknown samples. The variation observed in the composition of this material precipitated with the samples was approximately 0.2‰ (n = 8). Two additional materials, an inorganic apatite (HAP), and an in-house matrix-matching bioapatite enamel standard (BES) were also repeatedly used within the batches. The values obtained for these materials are 17.1 ± 0.3‰ (n = 6) and 17.0 ± 0.2‰ (n = 9), respectively, suggesting an overall good reproducibility of the precipitation method employed in this study.
As mentioned, because of the lack of an internationally recognised standard for phosphate oxygen isotope methods and commonly agreed protocols for precipitation and calibration, other laboratories might use different ways of correcting their δ 18 O p values. This makes comparisons of different datasets problematical. By using different calibration methods (1-point vs multiple point calibrations), for example, we can predict that the average possible shift between data calibrated in different ways can be between 0.3‰ and 0.5‰, but can reach 1‰, with Scientific RepoRts | 6:34986 | DOI: 10.1038/srep34986 the highest impact on the extreme values (i.e. those that have the most enriched or the most depleted δ 18 O compositions). For this reason, only data from the Beaker People Project were used to create the geostatistical model presented in this paper.