Approaching mercury distribution in burial environment using PLS-R modelling

Mercury environmental cycle and toxicology have been widely researched. Given the long history of mercury pollution, researching mercury trends in the past can help to understand its behaviour in the present. Archaeological skeletons have been found to be useful sources of information regarding mercury loads in the past. In our study we applied a soil multi-sampling approach in two burials dated to the 5th to 6th centuries AD. PLRS modelling was used to elucidate the factors controlling mercury distribution. The model explains 72% of mercury variance and suggests that mercury accumulation in the burial soils is the result of complex interactions. The decomposition of the bodies not only was the primary source of mercury to the soil but also responsible for the pedogenetic transformation of the sediments and the formation of soil components with the ability to retain mercury. The amount of soft tissues and bone mass also resulted in differences between burials, indicating that the skeletons were a primary/secondary source of mercury to the soil (i.e. temporary sink). Within burial variability seems to depend on the proximity of the soil to the thoracic area, where the main mercury target organs were located. We also conclude that, in coarse textured soils, as the ones studied in this investigation, the finer fraction (i.e. silt + clay) should be analysed, as it is the most reactive and the one with the higher potential to provide information on metal cycling and incipient soil processes. Finally, our study stresses the need to characterise the burial soil environment in order to fully understand the role of the interactions between soil and skeleton in mercury cycling in burial contexts.


Results
The fine earth fraction samples showed very low Hg concentrations (median = 2.6 ng g −1 , IQR = 1.7 ng g −1 , relative error average = 47%). Mercury concentrations in the silt + clay fraction varied from 14 to 24 ng g −1 for L01, 12 to 39 ng g −1 for L06 and 6 to 33 ng g −1 for L07 (Table 1). Thus, concentrations in the finer soil fractions (i.e. silt + clay) are up to tenfold the average of the fine earth mercury concentrations. Silt + clay mercury concentrations in samples outside the burials ranged between 6 and 23 ng g −1 and those inside the burials between 12 and 39 ng g −1 . www.nature.com/scientificreports/ Given the low concentrations in the fine earth, the PLSR analysis was performed on the silt + clay fraction, despite it represents a small part (0.26-12%) of the bulk soil. A summary of the values for the elements in the predictors' matrix can be found in Table 2. A PLSR model with three latent variables (LV) accounted for 72% of Hg variance ( Fig. 1 SM), corresponding 41% to LV1, 14% to LV2 and 17% to LV3 ( Table 3). The model regression coefficients were 0.32 β1, 0.31 β2, 0.73 β3 (Table 3). Figure 1 shows the relative weight of the LVs on the prediction of each sample. In general terms LV1 > LV2 > LV3, being the weight of LV3 relatively low for most of the samples. LV1 is characterised by positive loadings of Cu, Zn, and N and negative loading of Ca, Sr, and module (Table 4). LV2 shows positive loadings for P, C, Ca, and S and negative loadings for Ti and Fe (Table 4). LV3 is characterised by positive loadings for U, and S and negative loadings for P, and Mn (Table 4). Positive LV1 scores, were found for samples collected inside the burials, while negative scores were found for samples collected outside the burials (Fig. 2). LV2 positive scores mainly correspond to samples from burial T5, except for the most distal samples of the longitudinal transect of L07 (Fig. 2). LV3 scores show a more irregular pattern, showing positive scores on samples located in the thoracic area of the buried individuals (Fig. 2). The model residuals are low, with somewhat higher overestimation in the proximal samples of the L06 individual of T5 (Fig. 3).

Discussion
The low mercury concentrations of the fine earth samples are in agreement with previous studies in non-polluted soils [43][44][45][46] and are consistent with the predominantly sandy grain size of the samples and the absence of mercury mineral phases. The silt + clay fraction, on the other hand, is usually enriched in organic matter and metals 47 .
The PLSR model obtained for the silt + clay fraction provided three main latent vectors, i.e. underlying processes, that account for a significant proportion (> 72%) of the variation in mercury concentrations in the two burials. LV1, which explains the largest amount of mercury variability (Table 3 and Fig. 1), accounts for differences between inside and outside the burials (Fig. 2). LV2 can be related to differences between the two burials ( Fig. 2). LV3 seems to be due to within burial variability related to the proximity of the soil samples to areas where larger body mass (i.e. thoracic) was present (Fig. 2). Residuals showed no-clear pattern (Fig. 3) and can be linked to micro-scale soil variability.
The LV1 variables with positive loadings (Cu, Zn, N) are related to the abundance of organic matter in the silt + clay fraction. Copper and Zn are known to bind to the organic matter (OM) and clay 48 and N can be related to the total amount of OM 49,50 . Indicators with negative loadings include elements related to carbonates content (Ca and Sr) that can be linked to the presence of mollusc shells remains-note that the area is by the sea-and  www.nature.com/scientificreports/ location within the burial (module). The positive regression coefficient for this LV indicates that mercury content increases with organic matter content, which is clearly higher inside the burials (Fig. 4C). Many investigations have shown that mercury in soils and sediments is correlated to organic matter content 14,51 . This can be seen as a sort of fingerprint of the individuals that also affected the mercury content (Fig. 5). The buried bodies most probably were both the source of the organic matter (together with the wood coffin) and the source of the mercury, which has accumulated in the finest and most reactive soil fractions (i.e. OM and clay minerals). The original substrate is a dune dominated by silicic sand with biogenic carbonates that lacks the capacity to retain mercury. This is consistent with the low Hg concentrations in the samples outside the burial environment (both in the fine earth: median = 1.6, IQR = 0.8 ng g −1 ; and the silt + clay fraction: median = 12.1, IQR = 3.6 ng g −1 ). The LV2 indicators with positive loadings are P, C, and Ca, all constituents of the bone mineral component (i.e. hydroxyapatite), and S, which is a proxy of body decomposition. Among the indicators with negative loadings, Ti is indicative of soil fine fraction (i.e. clay) 52 . Iron is representative of major soil components (iron oxides) but also increases in the clay fraction as secondary oxy-hydroxides. LV2 scores are systematically positive in the double burial, T5, with the exception already mentioned before (distal samples in the L07 transect). This means that the fine soil fraction has higher concentrations of elements typical of the bone mineral matrix and S containing compounds than the same fraction in T1 burial. This difference most probably reflects the fact that T5 contains two individuals (an adolescent and a mature adult)-with a priori higher bone and soft tissue mass    www.nature.com/scientificreports/ interacting with the soil components-and T1 contains only one individual (a single burial). The weight of LV2 in the samples is systematically positive in L06 individual, while it is negative in L01 and shows an irregular pattern in L07 (Fig. 1). The regression coefficient for LV2 is positive meaning that mercury content is higher in T5 burial compared to T1. The link between mercury and the other elements can be chemical bone alteration and sulphur compounds availability during early body decomposition. Sulphur-containing amino acids (i.e. cysteine, cystine, and methionine) are subjected to desulphydralation by gastrointestinal microorganisms (i.e. sulphate-reducing bacteria) during body decomposition, producing the release of hydrogen sulphide gas, sulphides, ammonia, thiols, and pyruvic acids 49 . Given that during the main phase of tissues breakdown the tissue mass keeps anoxic 49 and Hg tends to precipitate as its sulphide form 51 , HgS compounds can be expected in these environments. 37 also hypothesised about this option, but provided no explanation for it. HgS compounds have extremely low solubility 51 , which will prevent from lixiviation, contributing to Hg retention in the burial environment. Thus, soft tissues decomposition may play a role in Hg release, as primary source (Fig. 5), but they are also key in providing compounds involved in Hg retention and accumulation in the burial environment.
Bone tissue acts as a sink for some elements during body decomposition 21 and then as a source to the burial environment, altering the elements that are released and accumulated in the finer, more reactive, soil fractions. Indeed, bone is a potential primary and/or secondary source of mercury to the soil of the burial; while soft tissues are the primary source both to soil and bone (Fig. 5). A fact that is in agreement with the recent finding that mercury can be stabilised in bioapatite as (Hg) 3 (PO 4 ) 2 27 . Furthermore, it could explain why mercury concentration in archaeological bones is higher than the concentration found in actual bones from autopsies and biopsies 53,54 -although lower contents in recent bone can also be due to lower exposure as today environmental loads are much more controlled. The higher the bone mass the higher the transfer of mercury to the soil. 42 also suggest that mercury distribution in the soil could be affected by burial context when graves are placed closed by as interactions between them cannot be excluded. In our case, T1 and T5 burials were far apart (~ 4 to 5 m) and with no other graves over or close to them. Another factor that should also be considered is that the overall mercury content of the individual in T1 compared to that of the individuals of T5 may have been different depending on environmental exposure during life. In a previous investigation, bone mercury content for individuals of the same period as those studied here varied between 4.6 and 299.9 ng g −113 .
The LV3 indicators with positive loadings are U and S. Uranium can be related to the geology of the area and has a high affinity for bone phosphates 55,56 , while S is a proxy for early stages of body decomposition, as previously mentioned in LV2 interpretation. Indicators with negative loadings are P, an indicator of bone alteration 55 , and Mn, a proxy for redox conditions 57 . The majority of samples with positive LV3 scores corresponds to those www.nature.com/scientificreports/ collected in the thoracic area of the three skeletons. The weight of LV3 tends to positive in samples of burial T5 while it does not show the same pattern in burial T1 (Fig. 1). This can be reflecting difference between burials. As already commented, the most important mercury target organs (kidneys, liver) are located close to the thoracic area 5 , which also contains the higher body mass-in comparison for example with the limbs or the head. Thus, differences in mercury content between soil close to thoracic area and soil close to the rest of the body would be expected. Furthermore, proteins forming the kidneys and liver start its putrefactive changes at early stages of body decomposition, while epidermis and muscles are more resistant to breakdown 49 . It is also likely that the higher content of Hg in this area is related with the higher availability during early body decomposition of both mercury and S compounds. Since the process takes place in the earliest stages of body decomposition the nonputrefied tissues would produce a confined environment and decrease mercury leaching. The negative loadings of P in this LV may also support this interpretation, as the soft tissue/bone ratio should also be higher in this body area and thus with lower net transfer of P to the soil. We speculate that "thoracic area" effect may be even more pronounced in graves with high mercury levels (i.e. periods with higher mercury levels of pollution and/ or individuals which suffered an acute intoxication, among other causes).
As already commented, the model residuals showed no clear pattern (Fig. 3) and can reflect soil microscale heterogeneity.
Diagenesis is of major concern when trace elements are studied in archaeological bones, as there is the potential for the soil to be a source of elements for bones 30,55 . But works which assessed mercury diagenetic incorporation into bone in archaeological burial environment 22,23,31,41,42 conclude that mercury bone levels were no related with diagenesis, unless the soils were naturally enriched in mercury. In a previous investigation in the same archaeological site we already reported very low mercury contents in soil (a pedo-sedimentary sequence produce an average of 0.9 ± 0.7 ng g −1 ) in comparison to bones of individuals of the same period as those of the present study (20.6 ± 23.5 ng g −1 ) 13 . Our PLSR model also supports the idea that the soil was not a significant source of Hg for buried bones. The situation seems to be rather the opposite; the buried individuals were most probably the source of the mercury accumulated in the soil (Fig. 5). In the particular case of our study, the decomposition of www.nature.com/scientificreports/ the buried bodies seems to have triggered a number of processes-developing a positive feedback mechanismthat lead to the rapid evolution of the parent soil material (i.e. sand dune), contributing to mineral weathering (i.e. decarbonation), clay formation and organic matter enrichment, creating the necessary conditions for the retention of the mercury released during the decomposition of the target organs (i.e. liver and kidney).

Conclusions
In the present research we approached mercury distribution in burial environment through PLSR modelling of elemental composition of soil/sediment samples. The fitted model explains 72% of mercury variability. We have obtained three main conclusions. i) Our data confirm human bodies (soft tissues and bones) as sources of mercury in burial environments, as well as sources of organic matter which also contributes to mercury fixation. Soil seems to have not been a source of mercury for the skeletons, since mercury levels in samples surrounding the burials were very low when compared with those inside burials. In addition, ii) burial context seems to affect mercury distribution. We found differences that could be explained by the number of individuals buried together, as the body mass (i.e.: soft tissues and bones) is related to mercury content in the burial soil. iii) The thoracic area seems to produce a further local enrichment in Hg, most probably due to its release from target organs and its retention in the reactive soil fractions (organic matter and clay). The obtained findings remark the importance of analysing soil/sediment associated to the burials in order to provide complementary information that can be key to understand some processes. Regarding the analytical procedure, we much encourage to analyse the finer (silt + clay) fractions when researching sandy soil/sediments, because this fraction is the most reactive and has better potential to contain information related to the body/soil interactions. In addition, the study of individuals/sites with low mercury levels can produce high quality results and complement those made in high polluted environment.
According to our results, soil/sediments from cemeteries can be significantly enriched in mercury compare with nearby ones. A fact that should be taken into account when describing soil properties or when those areas are having a new agrarian use. In addition, to assess the mercury content of an archaeological soil from a necropolis the proximity to the thoracic area of the skeletons can bias the obtained results.
Further research in mercury on soil/sediments associated to skeletons could help to achieve a better understanding of how mercury is distributed in these particular environments. Other questions such as how different burial contexts can affect mercury distribution in graves and how soft tissues affect it, remain open. Sulphur compounds, possibly related to body´s amino acids like methionine and cysteine, seem to contribute to Hg fixation-to analytically find these HgS compounds would be a necessary step in future research. HgS compounds are highly stable, limiting mercury release to the environment-even when the total amount is small. Therefore, to identify these compounds can help to prevent mercury release in individuals/cemeteries with high mercury concentrations.

Material and methods
Soil/sediment sampling. Soil and sediment samples of two post-Roman (AD 5th century) burials (T1 and T5) from A Lanzada site (Sanxenxo, Pontevedra, NW Spain, 42° 25′ 46″N 8° 52′ 25″W) (Fig. 4) have been analysed. A Lanzada is a well-known archaeological site with a widely studied necropolis (see among others 13,20,55,[58][59][60][61][62][63]. However, there is only an unpublished work about soil/sediment geochemistry 64 . Geologically, A Lanzada site is a dune that was subject to anthropic alteration-at least since the BC 2nd century-that reduced its dimensions 65 . The site was discovered in the 18th century and was mainly excavated during the 1960s and 1970s. Most of the effort during the earlier excavations was done in the necropolis. It is placed in the East side of the site and has two funerary areas-Roman and post-Roman 62 . Neither a pedological study nor sediment-soil sampling was developed during those campaigns and information about the soil characteristics where bodies were inhumated was based on photographs and colour changes in the bones (see 55 ). Skeletons were found in both dune sands and earth (i.e. acidic soils) independently of the funerary area 55 . Little effort was performed in the settlement area placed on the West side of the site, but recent campaigns in 2010 and 2016-2017 were mainly developed here.
The last excavation campaign took place in 2016-2017 and comprised the settlement area with houses occupied from the Late Bronze Age (~18 century BC) to the Late Iron Age-Early Roman Period (~ 1st century BC-AD). On the East margin of the excavation area a more monumental structure was discovered and dated to post-Roman times according to the archaeological material (~ 6th to 7th century AD). It was preliminary described as a church due to the monumental size of the walls 65 . Two coffin burials were discovered close to this structure and belonging to the same archaeological layer-named as T1 and T5. They were West-East oriented and radiocarbon dating indicates the individuals died between the 5th to the 6th centuries cal. AD. Both burials seem not to be related to the funerary areas previously found. T1 is a single burial containing a skeleton from an elder (> 60 years) female (L01). T5 is a multiple coetaneous burial containing two male skeletons, an adolescence (15-19 years) placed West-East (L06) and a mature adult (40-50 years) placed East-West (L07). All were in supine position with stretched arms and legs. Bodies were deposited presumably inside wood boxes (i.e. coffins) since nails were found surrounding the skeletons. The three skeletons were well-preserved, presenting the following abrasion degrees (according to 66 ): L01 degree 3, L06 degree 2, L07 degree 2. Despite L06 and L07 low superficial abrasion degree and general good preservation, these skeletons have intense physical (pressure) and chemical alteration in localised areas (L06: face and feet, L07: long bones epiphyses), which results in the loss of some of the bones/bone areas 67 .
Both burials were excavated in the dune. A clear colour pattern of rectangular shape was evident in both burials, marking the area delimited by the wood coffins (the nails were found at the edges). The colour is most probably related with increased organic matter content (Fig. 4C) www.nature.com/scientificreports/ and one transverse transect on each individual plus some additional samples taken from the inside the crania, and on L001 coxal area (Fig. 4). All samples have been taken on the field except for the additional ones, which were collected from soil adhered to the bones once in the laboratory. In total 46 soil/sediment samples were collected and analysed. Given that sample position is susceptible to affect mercury concentration, it was recorded using the longitudinal axis of each skeleton (Fig. 4). Sample coordinates were refereed to this reference axis and module was used as predictive variable.
Elemental composition analysis. Mercury content was analysed in two fractions, fine earth (< 2 mm, FE) and silt + clay (< 0.05 mm, SC). The samples were previously dried (25 °C) and finely milled (< 0.5 mm). Mercury concentrations were determined using a DMA-80 (Mileston) hosted at the laboratory of the Ecotoxicoloxía e Ecofisioloxía Vexetal research group (Departamento de Bioloxía Funcional, USC) following the protocol described in 13 . BCR277R (128 ± 17 ng g −1 ) was used as certified reference material. The quantification limit was 0.51 ng g −1 . Mean recovery for BCR27R was 91% (116 ± 17 ng g −1 ). Quality control has been done in 13 samples, reporting a relative error average of 6.37% for silt + clay. Elemental composition analysis was carried out on dried (25 °C) silt + clay samples using two energy dispersive X-ray fluorescence (EMMA-XRF) analysers 68,69 hosted at the RIAIDT facility of the Universidade de Santiago de Compostela (USC). This non-destructive technique provided concentration measurements for 10 elements (P, S, Ca, Ti, Mn, Fe, Cu, Zn, Sr, and U; see Table SM 1 for detection limits). Carbon and N analyses were done using a LECO (TruSpec CHNS) hosted at the Elemental Analysis Service of the RIAIDT in Campus Terra (USC).
Nitrogen was selected as a direct proxy for organic matter content (OM); Cu and Zn as indicators of metal enrichment in soil fine fractions bounded to OM and clay; Ca and Sr as proxies of biogenic carbonates, due to the presence of mollusc shells. The association of Ca, P and C is an indicator of the presence of bone, as these elements are the major constituents of bone mineral component (i.e. hydroxyapatite). Sulphur is used as proxy for OM forming strong bindings with Hg; Fe a major soil component; Ti is a proxy for fine grain fractions, and Mn is used as a proxy for red-ox conditions. Finally, U is characteristic of the geology of the area and presents high affinity for phosphate groups (i.e. bone mineral matrix).
Statistical methods: partial least squares regression modelling. To understand mercury distribution in the soil we used Partial Least Squares Regression (PLSR) 70 . The elemental composition data set (elemental concentrations in the silt + clay fraction; see explanation in the results section) was used as predictors while Hg concentration was set as response variable. Zero imputation has been performed for those cases where the element concentration was below the detection limit following the robust method proposed by 71 . Due to the compositional nature of our data the data set was transformed to centered log ratios (clr) 72 and scaled 70 . The main purpose of the PLSR modelling was the identification of the latent variables (i.e. processes) that affect Hg distribution and their weight in the soil samples. The model was performed using leave one out cross-validation without splitting the data in training and test sets to avoid the loss of relevant information and its explanatory power.
See SM for: PLSR data set, ŷ values, T matrix and W* matrix.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.