Structural equation modelling of mercury intra-skeletal variability on archaeological human remains

Structural equation modelling of mercury intra-skeletal variability on archaeological human remains Noemi Álvarez-Fernández a,b, ⁎, Antonio Martínez Cortizas a,c , Olalla López-Costas d,e,f a CRETUS, EcoPast (GI-1553), Facultade de Bioloxía, Universidade de Santiago de Compostela, 16782, Spain b Boscalia Technologies S.L., Spain c Bolin Centre for Climate Research, Stockholm University, Stockholm SE-10691, Sweden d EcoPast (GI-1553), CRETUS, Area of Archaeology, Department of History, Universidade de Santiago de Compostela, 15782, Spain e Archaeological Research Laboratory, Stockholm University, Wallenberglaboratoriet, SE-10691, Sweden f Laboratorio de Antropología Física, Facultad de Medicina, Universidad de Granada, 18012, Spain HIGHLIGHTS GRAPHICAL ABSTRACT •Skeletons have a dual role (sink and source) on Hg dynamics in graves. •Skeletal Hg variability seems to be affected by ante-mortem status of bone structure. •Bone crystallinity controls the role of bone components on Hg retention. •Soil plays a minor roleon Hg content in archaeological skeletons. ABSTRACTARTICLE INFO Editor: Mae Sexauer Gustin Keywords: Hg PLS-SEM Skeleton Soil/sediments Osteoarchaeology Diagenesis Archaeological burial environments are useful archives to investigate the long-term trends and the behaviour of mercury. In order to understand the relationship between mercury, skeletons and soil, we applied Partial Least Squares - Structural Equation Modelling (PLS-SEM) to a detailed, multisampling (n = 73 bone samples +37 soil samples) design of two archaeological graves dating to the 6th to 7th centuries CE (A Lanzada site, NW Spain). Mercury content was assessed using a DMA-80, and data about bone structure and the grave soil/sediments were obtained using FTIRATR spectroscopy. The theoretical model is supported by proxies of bone structure, grave soil/sediments, and location of the bone within the skeleton. The general model explained 61 % of mercury variance. Additionally, Partial Least Square –Prediction Oriented Segmentation (PLS-POS) was also used to check for segmentation in the dataset. POS revealed two group of samples depending on the bone phase (hydroxyapatite or collagen) controlling the Hg content, and the corresponding models explained 86 % and 76 % of Hg variance, respectively. The results suggestthat mercury behaviour in the graves is complex, and that mercury concentrations were influenced by i) the ante-mortem status of the bone matrix, related to the weight of each bone phase; ii) post-mortem evolution of bone crystallinity, where bone loses mercury with increasing alteration; and iii) the proximity of the skeletal pieces to mercury target organs, as Science of the Total Environment 851 (2022) 158015 ⁎Corresponding author at: CRETUS, EcoPast (GI-1553), Universidade de Santiago de Compostela, 16782, Spain. E-mail address: [email protected] (N. Álvarez-Fernández). http://dx.doi.org/10.1016/j.scitotenv.2022.158015 Received 31 May 2022; Received in revised form 2 August 2022; Accepted 9 August 2022 Available online 13 August 2022 0048-9697/©2022 TheAuthors. Published by Elsevier B.V. Thisis an openaccess articleunder theCC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4. 0/). Contents lists available at ScienceDirect Science of the Total Environment journal homepage: www.elsevier.com/locate/scitotenv decomposition and collapse of the thoracic and abdominal soft tissues causes a secondary mercury enrichment in bones from the body trunk during early post-mortem. Skeletons provide a source of mercury to the soil whereas soil/ sediments contribute little to skeletal mercury content. 1. Introduction The cycle and toxicology of mercury are widely researched due to its toxicity. Mercury is a global public health concern according to the World Health Organization (WHO, 2020) and, in general terms, it has been regarded as the most toxic non-radioactive element (Pushie et al., 2014). It is harmful even at very low doses, and has no known biological function. Furthermore, mercury is a worldwide distributed pollutant. There are both anthropogenic and non-anthropogenic sources to the environment, and once it is bioavailable mercury bioaccumulates and biomagnifies in food chains and is very persistent in ecosystems (Evers, 2018;Morel et al., 1998;Tang et al., 2020). Mercury has different chemical fractions that occur naturally in the environment; including elemental mercury (Hg 0 ), inorganic mercurous (Hg + ) and mercuric (Hg 2+ ) salts, as well as organic compounds such as methyl- and ethyl‑mercury (Berlin et al., 2015). The primary toxic effects of this metal are caused by its capacity to bind to sulfhydryl groups, and, to a lesser extent, to hydroxyl, carboxyl, and phosphorylgroups. This ability to bind such compounds affects non-proteins containing thiol groups, reduced glutathione and proteins containing cysteine, in addition to its interactions with ionic channels, transporters, and enzymes (Silva-Filho et al., 2021;Tchounwou et al., 2003). Mercury's different chemical forms do not necessarily share absorption paths or behaviour once they are inside the organism, which makes it toxicologically complex (Berlin et al., 2015; Clarkson, 1997). Furthermore, toxicity is also dependent on the dose, time, and way of exposure (Abass et al., 2018). This variability leads to a wide range of target organs and clinical symptomatology. At low-dose and chronic exposure (i.e. environmental) kidneys and liver are the main target organs (García et al., 2001). As a result of efforts to eliminate the direct sources of mercury pollution, which in some cases have led to acute poisoning cases, most studies on mercury toxicology focused upon the effects of high-dose exposures (see Álvarez-Fernández et al., 2020); despite cases are predominantly of lowdose (Holmes et al., 2009). Chronic exposure at low doses is primarily related to anthropogenic emissions and the resulting environmental pollution. In these circumstances (low dose and chronic), mercury has important consequences for health, such as chronic intoxication and increased risk of developmental disorders in children (Budnik and Casteleyn, 2019;Ha et al., 2017). Understanding the mercury cycle is a subject of intensive research recently (Gustin et al., 2020 and papers contained therein). This cycle is altered by the long residence time of mercury in natural systems (Outridge et al., 2018;Streets et al., 2011). Currently, 60 % of the annual emissions to the atmosphere come from mercury previously deposited on soil and water (Li et al., 2022). Mercury in soils is mainly bond to organic matter, a temporary sink until it is re-mobilised (Gabriel and Williamson, 2004;Gębka et al., 2020;Skyllberg et al., 2003). Cemeteries are a sources of soil mercury, especially when the inhumation of the body was part of the funerary ritual (Amuno, 2013;Jonker and Olivier, 2012;Mohammed and Abudeif, 2020;Prestes da Silva et al., 2020;Spongberg and Becks, 2000;Uslu et al., 2009;WHO, 1998). Human bodies act as temporary sinks during life as they incorporate mercury from the environment by inhalation, ingestion, etc. (Liu et al., 2011). Once bodies are buried, they act as a source of mercury –among other elements –to the soil, and leave afingerprint on it (García-López et al., 2022;Graf, 1986;Sobocká, 2004). When a large number of bodies are buried in the same place (i.e. longused cemeteries), the transfer of mercury is significant and could have an adverse impact on surrounding ecosystems. Natural archives show that humans have been impacting the mercury cycle since at least 3250 BCE (e.g., in the South Iberian Peninsula) (Leblanc et al., 2000). In NW Spain atmospheric mercury pollution has been dated back to 500 BCE (Martınez-Cortizas et al., 1999), and can be related to mercury mining in the Iberian Peninsula during the Iron Age. Anthropogenic emissions markedly increased to peak during the Roman Empire due to the extensive mining and metallurgy to fulfil economic demands. Pollution decreased with the fall of the Western Roman Empire. Álvarez-Fernández et al. (2020) and López-Costas et al. (2020) assessed the impact of these fluctuations on mercury levels in past-populations through their osteoarchaeological remains, showing similar trends to those provided by local reconstructions of atmospheric pollution from peatland archives. Therefore, mercury contamination in skeletons can complement the data obtained from commonly used natural archives and provide more specific information about the impact of certain pollutants, like mercury, on biota (Cooke et al., 2020;Gustin et al., 2020). However, the preserved archaeological remains are mainly bones and teeth (i.e. skeletons), which are not the main target organs for mercury and, in some cases, have been buried for a long time (for a review of Hg studies on archaeological populations see Álvarez-Fernández et al. (2020)). Diagenesis is a major issue when trace elemental composition is studied in osteoarchaeological records (Hedges, 2002). For the graves studied in this work, we addressed the factors controlling mercury content in soil/ sediments in a previous investigation (Álvarez-Fernández et al., 2021). This study showed that the buried individuals were the only source of mercury, and results agreed with previous work that found no evidence of diagenetic incorporation of mercury from soil to bone (Emslie et al., 2015;Kepa et al., 2012;Rasmussen et al., 2008, 2013b;Walser et al., 2019;Yamada et al., 1995). However, the potential impact that past human demography and funerary rituals may have had on the soil current state remains to be investigated. To do so, it is important to address the relationship between mercury and bone tissue, both ante- and post-mortem, to be able to understand the role of skeletons on mercury dynamics in the surrounding soil. In this investigation, two post-Roman-Early-Medieval burials from the archaeological site of A Lanzada (NW Spain) were selected and a detailed soil and bone sampling scheme wasused to assess the variability of mercury concentrations in the skeletons and the soil. Bone samples of the buried individuals were taken as close as possible to the location of the soil samples. The weight of potential factors determining the intra-skeletal variability in mercury concentrations was assessed using Partial Least Squares - Structural Equation Modelling (PLS-SEM), with the objective of elucidating: i) the role of bone structure (ante- and post-mortem)inmercury content; ii) the role of bone components on mercury retention; iii) whether the bones were a source and/or sink of mercury in the burial environment; and iv) the role of skeletons on the mercury cycle. 2. Material and methods 2.1. Bone and soil samples Seventy-three bulk bone and 37 soil/sediment (silt+clay fraction) samples associated with three skeletons found in two post-Roman (AD 5th century) burials (T1 and T5) from A Lanzada site (Sanxenxo Potenvedra, NW Spain, 42° 25′46″N8°52′25″W) were analysed (Fig. 1). A Lanzada is an archaeological site with several occupational phases. The most relevant ones to this study are a long-used habitational area (Bronze to post-Roman times), to the West, and a necropolis with two funerary areas (Roman and post-Roman; excavated during the 1960s and 1970s) to the East (López-Costas, 2015;Rodríguez Martínez, 2017). The graves, T1 and T5, were included in this study, were discovered on the West side of the site during the last excavation campaign (2016–2017), when a detailed anthropological and pedological study was N. Álvarez-Fernández et al. Science of the Total Environment 851 (2022) 158015 2 done (including the collection of soil samples from the burials). They were located closed to each other and, according to the archaeological material (~6th to 7th century CE), dated to the same archaeological layer of a monumental structure that was thought to be a church (Rodríguez Martínez, 2017). The burials were oriented West-East and excavated from dune sands. The grain size of the burials soil/sediments is dominated by sand (85 % in the burials' soil and 95 % in the soil outside the burials; García- López et al., 2022). The burial soils showed an enrichment in P and a decrease in alkalinity (from pH ~9 to 8), compared to the soil located from outside the burials. They also showed a secondary enrichment in silt+clay fractions and organic matter (Álvarez-Fernández et al., 2021), soil components able to retain mercury. T1 and T5 are not directly related to the funerary areas of the necropolis previously found, so the large number of bodies there was not expected to greatly influence the geochemistry of the two burials. Radiocarbon dating indicated that the individuals died sometime between the 5th to the 6th centuries cal. AD. Bodies were probably deposited in wood boxes (i.e., coffins) since a clear colour pattern of rectangular shape was evident in both burials (Fig. 1) and nails were found at the edges. All the skeletons were in supine position with stretched arms and legs (hands placed at both sides of the body or over the pelvic area; see Fig. 1). T1 is a single burial containing the skeleton of an elderly (>60 years) female (L01) of short stature (152 cm) when compared to Iberian archaeological populations from Roman and Medieval times (López Costas, 2012;López Costas et al., 2017). T5 is a multiple coetaneous burial containing two male skeletons, an adolescent (15–19 years) placedin West-East orientation (L06) and a mature adult (40–50 years) placed in East-West orientation (L07). Both had above-average statures (175 cm and 170 cm, respectively) when compared to archaeological populations from the same period (López Costas, 2012;López Costas et al., 2017). Sex was estimated following established criteria on the innominate and cranial bones (see a summary in Buikstra and Ubelaker (1994)). Age was estimated using standard identification criteria for innominate bone (auricular surface and pubic symphysis morphology), fourthrib and epiphyseal fusion (for a summary see Buikstra and Ubelaker (1994)), and Iberian methods for growth and maturity of postcranial bones (López-Costas et al., 2012; Rissech et al., 2013). Statures were calculated using the humerus and femur maximum length (De Mendonça, 2000). All bone remains were well-preserved (Table SM 1), presenting the following abrasion degrees (average for the whole preserved bony parts, according to Brickley and McKinley, 2004): L01 degree 3, L06 degree 2, L07 degree 2. L06 and L07 have intense physical (pressure) and chemical alteration in localised areas (L06: face and feet, L07: long bones epiphyses), that resulted in the loss of some of the bones/bone areas despite their low degree of superficial abrasion and being well preserved (López Costas et al., 2017). Bulk bone (cortical) was micro-sampled using a dentist's drill, after the superficial layer was removed to avoid contamination. Due to the importance of preserving archaeological skeletons, a small amount was sampled (~30 to 60 mg). Four types of bone were analysed when possible: i) type 0a: vertebrae (spine) and ilium; ii) type 1a: ribs; iii) type 2a: long bones; and iv) type 3a: crania (Table SM 1). This classification is modified from previous studies investigating intra-skeletal variability of mercury and other elements from the large necropolis area from A Lanzada site (Álvarez-Fernández et al., 2020;López-Costas et al., 2016). The sampling strategy followed a multisampling approach to cover most of the skeleton intra-variability (Fig. 2 and Table SM 1). For ethical reasons the number of samples in L01 was reduced due to the low bone density, related to severe osteopenia (possibly a case of senile osteoporosis). A total of 73 bulk powder bone samples were analysed: 11 for L01, 31 for L06 and 31 for L07. Since the present research is the continuation of previous work (Álvarez-Fernández et al., 2021), in which the distribution of mercury in the burial soil/sediment was assessed for these two graves, bulk bone sampling was intended to be in harmony with the soil/sediment sampling design, where samples were collected in two transects –longitudinal and transverse –along each skeleton. The soil/sediment samples used here and in the previous research are the same, excluding those from outside burials that were not included in this study. As soil/sediment samples were collected in the field, an exact match was not possible for all the bone samples, so the closest soil/sediment sample to each bone sample 50 m 5 cm 5 cm L001 L006 L007 T1 T5 A B C Fig. 1. (A) and (B) Aerial View of A Lanzada with approximated location of the graves; (A) modified from Google Earth 2020., (B) modified from Rodríguez Martínez, 2017. (C) Tombs. N. Álvarez-Fernández et al. Science of the Total Environment 851 (2022) 158015 3 was used. Median distance between bone and soil samples was <10 cm in all cases. Unfortunately, the other skeletons belonging to the large necropolis area were excavated >50 years ago without collecting soil samples. 2.2. Mercury analyses Bulk bone 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 Álvarez-Fernández et al. (2020). Two standard reference materials were used: estuarine sediment GBW07601a (670 ± 60 ng g −1 ) and bone meal SRM 1486 (2.3 ± 1.4 ng g −1 ). The quantification limit was 0.26 ng g −1 . Mean recovery for reference materials was 107 % GBW07601a (698 ± 50 ng g −1 ) and 104 % SRM 1486 (2.5 ± 1.2 ng g −1 ). Quality controls (replicated analyses every 10 samples) were done for 29 samples, reporting an average difference of 1.9 ng g −1 . 2.3. FTIR-ATR analyses Bone and soil samples were analysed by total attenuated reflectance Fourier-transform infrared spectroscopy (FTIR-ATR) using a spectrometer Agilent Cary 630 FTIR coupled with an ATR module, located at EcoPast laboratories (Facultade de Bioloxía, USC). Spectra were acquired in the mid-infrared region (MIR) 4000–400 cm −1 , by averaging 100 scans at a resolution of 4 cm −1 . The equipment was cleaned, and a background was collected before every measurement. Peak identification was done using the {andurinha} R package (Álvarez Fernández and Martínez Cortizas, 2020) following the second derivative sum spectrum method after Z score standardisation. Assignment of compounds related to vibrations is based on those reported in the literature (see references in Table SM 2), taking into account the limitations imposed on IR interpretation of complex samples (Coates, 2000;Larkin, 2017;Simonescu, 2012;Socrates, 2004). Several indices related to the bone structure were calculated: i) the infrared splitting factor (IRSF), which is related to hydroxyapatite crystallinity (Weiner and Bar-Yosef, 1990), as the sum of the ν 4 (PO 4 ) peaks intensities at 600 and 559 cm −1 divided by the intensity of the valley between them (587 cm −1 ); ii) the mineral maturity index (MMI), which represents the progressive transformation of poorly crystallised non-apatite domains into well crystallised apatite (Farlay et al., 2010), as the ratio 1019/ 1111 cm −1 of ν 3 (PO4); iii) the CO 3 /PO 4 (CP) ratio, which indicates the carbonate content of hydroxyapatite, dividing the ν 3 (CO 3 ) peak intensity at 1409 cm −1 by the ν 3 (PO4) peak intensity at 1019 cm −1 (Grunenwald et al., 2014); and iv) the collagen content estimated through the Amide I/ PO 4 ratio (AmP), dividing the peak intensity of ν1(Amide I)bandat 1638 cm −1 by the ν 3 (PO 4 ) band at 1019 cm −1 (Trueman et al., 2004a). 2.4. Statistics 2.4.1. Descriptive methods Statistical analyses and graphs were done using the software R (R Core Team, 2021). Normality was assessed using the Shapiro-Wilk test, as it is the most robust even when the sample size is relatively small (Yap and Sim, 2011). Given the non-parametric nature of the data, descriptive analysis was based on median and interquartile ranks (IQR). To detect significant differences among groups, Mann-Whitney-Wilcoxon was applied when comparing two groups, or Kruskal-Wallis combined with pairwise Mann-Whitney-Wilcoxon when more than two groups were involved. The p-values were considered significant when p <0.05. 2.4.2. Partial least squares structural equation modelling To assess mercury intra-skeletal distribution, Partial Least Squares Structural Equation Modelling (PLS-SEM) was applied, using the SmartPLS software (Ringle et al., 2015). This method enables the estimation of complex models following a causal-predictive approach without imposing distributional assumptions on the data matrix (Hair et al., 2019a;Sarstedt et al., 2017) while aiming to maximise the explained variance of the dependent latent constructs (Hair et al., 2019b). The inner and outer models are represented in Fig. 3. The latent constructs can be classified into three groups: i) those related to bone components, collagen and hydroxyapatite (HAP), to assess their role on the skeleton mercury content, and the degree of crystallinity of the mineral phase (HAPc), as it can affect the release of mercury from the bone tissue to the surrounding soil; ii) soil components Fig. 2. Sampling design. Sampled bones coloured in light purple with the sampling area highlighted with dots in dark purple. See further information in Table SM 1. N. Álvarez-Fernández et al. Science of the Total Environment 851 (2022) 158015 4 in the silt+clay fraction, among which mercury (Hg), primary silicates (psilicates) and clay were selected to assess the degree of interaction between bone and soil, as the bone was found to be a potential source of mercury to the soil; and iii) location, as it is susceptible to affect mercury content (Álvarez-Fernández et al., 2021). Measured variables included: i) for collagen the FTIR-ATR bands at 3280 (Amide A and –OH vibration), 1638 (Amide I) and 1545 cm −1 (Amide II) (Goormaghtigh et al., 2006;Martínez Cortizas and López- Costas, 2020;Trueman et al., 2004b); ii) for HAP bands at 870 (ν 2 (CO 3 )), 712 (ν 4 (CO 3 )), 690 (ν 4 (PO 4 ) and –OH), and 470 cm −1 (ν 2 (PO 4 )) (Baxter et al., 1966;Fleet, 2009;Rey et al., 2011;Zammel et al., 2021); iii) for HAPc the band at 1019 cm −1 characteristic of ν 3 (PO4) (Dal Sasso et al., 2016;Rey et al., 1991); iv) for soil Hg, concentrations obtained in a previous work on the same graves were used (Álvarez-Fernández et al., 2021); v) for p-silicates the bands at 606 (Si, Al tetrahedral deformation), 585 (Si tetrahedral breathing) and 531 cm −1 (Si –O–H vibrations) (McKeown, 2005;Pérez-Rodríguez et al., 2016); vi) for clay, bands at 3694 and 3621 cm −1 characteristic of kaolinite –OH stretching (Vaculíková et al., 2011); and vii) for location, the longitudinal axis (xaxis) of each skeleton setting as zero the pelvic area (approximately at S3 sacral vertebra). All measured variables were set on reflective mode. Given the compositional nature of the manifest variables for latent constructs bone Hg and soilHg, concentrations were previously transformed by natural logarithm. The data matrix (Table SM 3) was standardised and the path method was chosen for the weighting scheme (Hair et al., 2017). Twotailed bootstrapping was applied to test the statistical significance of path coefficients choosing the bias-corrected and accelerated (BCa) method for the confidence intervals (Hair et al., 2017). The significance level was set to 0.05. Partial Least Square Prediction-Oriented Segmentation (PLS-POS) was also used to identify segments (i.e., groups of samples) in the data set (Becker et al., 2013). 3. Results 3.1. Mercury variability Bone mercury concentrations varied from 2 to 4 ng g −1 for L01, from 1 to 39 ng g −1 forL06,andfrom1to18ngg −1 for L07 (Table 1). Differences in mercury between individuals were significant (p-value <0.01). L01 bones were the least enriched in mercury followed by L07, while L06 had the higher mercury concentrations. To evaluate mercury variability among bone types the four groups were considered. No significant differences were found between type 0a (vertebrae+ilium) and type 1a (ribs) and between type 2a (long bones) and type 3a (crania). Therefore, these groups were combined into two final groups (type1 = 0a + 1a: spine, ilium, ribs; and type2 = 2a + 3a: long bones, crania). From here on we will refer exclusively to these two bone types. Mercury concentrations Fig. 3. PLS-SEM inner model and outer model. HAP: hydroxyapatite; HAPc: hydroxyapatite crystallinity; p-silicates: primary silicates. Table 1 Summary of Hg concentrations in bone per individual and bone type (ng g −1 ). L01 is a female and L06 and L07 are males. Minimum Median Maximum IQR Individual L01 1.97 2.94 4.31 0.94 L06 1.19 9.21 38.58 11.69 L07 1.48 5.72 17.86 3.32 Bone type Type 1 2.15 8.92 38.58 8.35 Type 2 1.19 3.34 17.86 2.45 N. Álvarez-Fernández et al. Science of the Total Environment 851 (2022) 158015 5 were found to be statistically higher (p-value <0.01) in bone type 1 (varying from 2 to 39 ng g −1 ) than in bone type 2 (varying from 1 to 18 ng g −1 )(Table 1). 3.2. FTIR-ATR spectra Fig. SM 1 gives an overview of the full bone spectra in the range 4000–400 cm −1 . The average spectrum shows the most characteristic peaks, all the spectra enable to detect the regions with the highest variability among samples, while the second derivative spectrum locates peaks, including those with lower absorbance. The clay region, ~3700 to ~3610 cm −1 , showed high variability although average band intensity is low. The region from ~3600 to ~2800 cm −1 , with bands related to collagen vibrations (amide A and B, aliphatics and –OH groups), showed moderate variability and low absorbance. The Amide I collagen band, at ~1660 cm −1 , showed moderate absorbance but high variability. Absorptions in the region ~1450 to ~1300 cm −1 , corresponding to the ν 3 (CO 3 ) vibration modes, showed moderate intensity and variability. From ~1200 to ~950 cm −1 , the region with the characteristic ν 3 (PO 4 )andν 1 (PO 4 ) vibrational modes of phosphates, showed high to moderate intensities and low variability. Bands between ~900 and ~700 cm −1 ,ofν 2 (CO 3 ) and ν 4 (CO 3 ) vibrational modes of carbonates, showed moderate intensities and variability. The typical doublet of the ν 4 (PO 4 ) vibrational mode of phosphate, between ~690 and ~500 cm −1 , showed high intensities, but low variability. 3.3. Bone structure indices Bone mercury was slightly correlated with the infrared splitting factor (IRSF), mineral maturity index (MMI), carbon/phosphate ratio (CP), and amide I/phosphate ratio (AmP). The correlation with IRSF and MMI was negative with values around ~−0.36 (Table SM 4). Correlation with CP and AmP was positive, with values of 0.34 and 0.50 respectively. Variability was minimal when comparing individuals IRSF (L01: 4.3 ± 0.3, L06: 3.9 ± 0.9, L07: 4.1 ± 0.5), MMI (L01: 9.0 ± 0.9, L06: 8.2 ± 1.0, L07: 9.1 ± 2.0) and CP (L01: 0.07 ± 0.02, L06: 0.09 ± 0.06, L07: 0.09 ± 0.04). While for AmP significant differences were found between T1 (L01: −0.03 ± 0.00) and T5 (L06: 0.09 ± 0.06, L07: 0.09 ± 0.04) (Table 2). Similar results were obtained when comparing bone type, with low variability in IRSF (type 1: 4.0 ± 0.7, type 2: 4.2 ± 0.6), MMI (type 1: 8.5 ± 1.2, type 2: 9.3 ± 1.7), CP (type 1: 0.09 ± 0.04, type 2: 0.08 ± 0.03) and significant differences in AmP (varying from −0.03 to 0.05 in bone type 1 and from −0.03 to −0.00 in bone type 2). The bone type 2 sample L06.30 was identified as an outlier for AmP with a value of 0.07. 3.4. PLS-SEM model A PLS-SEM model (Fig. 3) with 8 latent constructs (i.e., inner model) and 16 indicator variables (i.e., outer model) was tested in the whole data set (G0). This general model (G0) accounted for 61 % of Hg variance. Additionally, PLS-POS analysis revealed the existence of two groups (Table SM 3). The same model (Fig. 3) used for G0 was applied to the two groups (G1 and G2). The model for G1 accounted for 86 %, and the model for G2 accounted for 76 % (Table SM 7 and Fig. SM 2) of the mercury variance. Interactions between LVs (i.e., paths)in G0accountfor 62 % of collagen variance, 59 % in G1 and 66 % in G2. The models also explain ~90 % of hydroxyapatite (HAP) variance (93 % G0, 95 % G1, 92 % G2). G0 also accounts for 38 % of soil Hg variability, G1 accounts for 17 % and G2 for 60 %. Table SM 8 shows the path coefficients for G0, G1 and G2 with the significance levels. The path coefficients from HAP to bone Hg (HAP➔bone Hg) are significant in the three models and high and positive in G0, moderate/high and negative in G1 and strong and positive in G2. The path coefficients hydroxyapatite crystallinity (HAPc)➔bone collagen and HAP were significant and negative in all the models; while the path coefficients HAPc➔bone Hg were high and negative for G1 and positive for G2, but not significant for G0. Bone collagen➔bone Hg was the only significant for G1, in which the path coefficient is moderate and positive. Location➔ bone Hg path coefficients were low and negative for the three models, being significant only for G0. Soil clay had a low and negative effect on bone Hg in the three models. Soil primary silicates (p-silicates) had a negative coefficientsin all the models, being moderate and significant only in G0 and G2. Soil clay showed very low and non-significant path coefficients on soil Hg in all the models, while soil primary silicates had a negative path coefficient in all of them, but significant only for G0 and G2. Finally, bone Hg has a low positive path coefficient on soil Hg in all the models, which were significant for G0 and G2. Fig. 4 shows the direct and indirect effects of all the model paths (Table SM 9). Soil p-silicates had a negative effect (direct and indirect) both on bone Hg and soil Hg in all models including: i) Location which has a negative indirect effect on soil Hg and a negative direct effect on bone Hg; ii) Bone collagen had a direct positive effect on bone Hg, but it was only significant for G1; iii) Hydroxyapatite crystallinity (HAPc) had a negative indirect effect on soil Hg and a negative indirect effect on HAP and bone collagen; iv) HAPc had a positive effect on bone Hg, being direct for G0 and G2 and indirect for G1. HAPc also had a negative effect on bone Hg, indirect in G0 and G1 and direct in G2; v) HAP had an indirect effect on soil Hg, being positive in G0 and G2 and negative, but non-significant for G1; vi) HAP had a direct effect on bone Hg, positive in G0 and G1 and negative in G1; vii) Bone Hg had a positive direct effect on soil Hg. The G0 residuals (Fig. SM 4) for L06 and L07 showed an underestimation of the samples in the thoracic area while overestimates the samples from long bones samples; no clear pattern was observed for L01 residuals. In G1 and G2 residuals were low and show no clear spatial pattern. Samples assigned to G1 and G2 did not respond to significant differences in Hg content, bone structure indices, or to differences between graves and among individuals, but a possible trend in the classification based on a visu of the skeletons was identified. Samples in group-1 were less ossified (thinner cortical) and had higher chemical alteration (surface abrasion) than those in group-2 (Table SM 1). Further information about the model can be found in the Supplementary Material, including block unidimensionality, outer model summary, inner model summary, inner model paths summary and total/direct/indirect effects. 4. Discussion 4.1. Mercury content Mercury content varied among individuals and decreased with age (L01 <L07 <L06). However, the small sample size (n = 3 individuals) does not allow us to discard a spurious relationship, as many other factors (e.g., exposure during life) could also be responsible for this pattern. In a previous study, the impact of mercury atmospheric concentrations in A Lanzada Roman (n = 43) and post-Roman (n = 33) population was assessed by analysing three types of bones (ribs, long bones, and crania). No relationship was found between mercury content and age in any of the two cohorts (Álvarez-Fernández et al., 2020), in line with results from recent research from biopsies (Babuśka-Roczniak et al., 2021;Zioła- Frankowska et al., 2017) and autopsies (Domingo et al., 2017;Yoo et al., 2002). Variation in mercury concentrations among the three analysed individuals (L01, L06, and L07) may be explained by individual differences Table 2 Summary of AmP concentration in bone per individual and bone type. L01 is a female and L06 and L07 are males. Minimum Median Maximum IQR Individual L01 −0.028 −0.025 −0.021 0.003 L06 −0.033 −0.016 0.071 0.025 L07 −0.030 −0.021 0.051 0.010 Bone type Type 1 −0.033 −0.015 0.052 0.022 Type 2 −0.031 −0.024 −0.002 0.004 N. Álvarez-Fernández et al. Science of the Total Environment 851 (2022) 158015 6 Fig. 4. PLS-SEM models: direct and indirect effects; significant paths are highlighted with (*). N. Álvarez-Fernández et al. Science of the Total Environment 851 (2022) 158015 7 in ante-mortem exposure and accumulation. Previous research at A Lanzada suggested that mercury exposure was mainly related to atmospheric pollution during both Roman and post-Roman times (Álvarez-Fernández et al., 2020;López-Costas et al., 2020). Exposure to larger doses of mercury was very unlikely as the median concentration was ~6 ng g −1 for in L01, L06, andL07and21±23ngg −1 for the post-Roman cohort from the necropolis. L06 and L07 were likely to share a similar environmental exposure, as they were buried in a double non-consecutive burial that indicates that both died at the same time or within a few days of each other. L01 could have lived a few decades earlier or later, and this individual could have had different environmental exposure. Differences in environmental exposure could explain in part the slight differences observed here, as well as those found in the individuals buried in the funerary areas of the same archaeological site (Álvarez-Fernández et al., 2020). In contrast, mercury accumulation appears to be related to individual factors such as daily life habits, occupational activities, and body mass index (Bjørklund et al., 2017;Zioła-Frankowska et al., 2017). So, differences in mercury content among the studied individuals may be explained by non-measured characteristics, such as body mass. If stature is used as a proxy for body mass (notwithstanding any limitations), the same trend as that observed in mercury content is found: L01 (152 cm) <L07 (170 cm) <L06 (175 cm). L06 and L07, both males, had higher stature, taller than the Iberian average for the period, thus may have had higher body mass than L01, a senile woman of shorter stature. L01 also suffered from advanced osteopenia, probably as a side effect of her age (>60 years-old) that could have affected mercury accumulation/release (Lanocha et al., 2013;Zioła-Frankowska et al., 2017). Old age could also have prevented her from participating in active occupational activities, and she may have had a more sedentary lifestyle compared to the two men. There is an ongoing debate regarding the association between mercury and osteopenia with contrasting results. Some studies suggest mercury affects bone metabolism (Suzuki et al., 2004; Tang et al., 2022), while others did not find any significant association between mercury exposure and this pathology (see a review in Jalili et al., 2020). There is only one study where mercury was found to disrupt calcium metabolism, in goldfish scales, and mercury dose was high (10 −7 Mmethylmercury exposure during 2, 4, and 8 days in Suzuki et al. (2004)). It is unlikely that mercury exposure was the cause of osteopenia in L01. Bone type showed significant differences in mercury content. Bone type 1 (spine+ilium+ribs) had higher mercury concentration than bone type 2 (long bones+crania). Rasmussen et al. (2017) found that ribs had a higher mercury content compared to long bones when individuals were treated with mercury remedies, and they proposed two causes: i) the bone turnover, which is quicker in short bones, and ii) the collapse of Hg-enriched soft tissues contained in the thorax. Both causes are possible and may coexist as ante- and post-mortem mechanisms, and therefore account for mercury variability between the two analysed bone types. Type 1 bones also have a priori higher turnover rates than type 2 bones, particularly long bones such as the femur (Hedges et al., 2007). Therefore, the higher Hg content in bone type 1 may be related to increased exposure to mercury in the most recent years before death, and possibly the main cause of variability. A large proportion of an adult femur is formed during growth spurts (Hedges et al., 2007), and it was produced when the individual was possibly less engaged in activities related to metal exposure (i.e., metallurgy). Emslie et al. (2019, 2015) observed that humeri had higher mercury content than femora and tibiae. They hypothesised that higher mercury content may be due to bone remodelling rates, relating themto biomechanical regulation as consequence of the use of heavy tools and other daily life activities. However, they did not explain why skeletal markers were more expressed in humeri than on lower limbs. Since they only analysed long bones (humeri, femora, and tibiae), the comparison with small bones from thoracic and abdominal area was not possible. An alternative hypothesis can be considered: the secondary mercury enrichment of bones from the body trunk occurred during early post-mortem, due to the decomposition and collapse of the thoracic and abdominal soft tissues (note that during the breakdown of soft tissues the body mass creates anoxic conditions (Janaway et al., 2009)). This finding is consistent with Álvarez-Fernández et al. (2021), who found that mercury content in the soil/sediments from these two burials studied here was related to proximity to the thoracic/abdominal area. The main target organs for mercury, when exposure is low and chronic, are placed in the abdomen and thorax (i.e.: kidneys and liver) and that lungs and digestive apparatus can also accumulate some mercury (Lech and Sadlik, 2004). Rasmussen et al. (2013a, 2013b) also measured mercury content in the soil close to the kidneys, liver and lung (n kidneys =3,n liver =4,n lungs = 4 individuals) revealing the same trend. For the whole population and archaeological/historical framework of A Lanzada, no evidence was found of direct use of mercury that caused an acute intoxication/high level exposure. On the contrary, it seems to have been low and chronic, through atmospheric pollution (Álvarez-Fernández et al., 2020); a pathway that is more likely to explain a higher proportion of the mercury variability than the turnover rate in the individuals studied here. Rasmussen et al. (2017), also suggested that differences in mercury content for individuals were due to treatment with mercury-containing remedies, with mercury concentrations between 80 and 78,727 ng g −1 in cortical bone. Several studies reported values in kidneys and liver from populations exposed to mercury. Johansen et al. (2007) found values of 1400 ng g −1 in kidneys, and 540 ng g −1 in liver (in Greenland populations with high fish consumption), and García et al. (2001) reported values of 250 ng g −1 in kidneys and 140 ng g −1 in the liver (from a Spanish population living near to an industrial area). Thus, the higher mercury content on the thoracic and abdominal region may be a proxy fingerprint of the soft tissues placed there, which include the main target organs for mercury when exposure is low and chronic (Holmes et al., 2009). The influence of the thoracic andabdominal area could be extendedto the long bones placed near them, especially the humerus. Depending on the buried body position, ulna and radius may be parallel to the legs rather than over the abdomen (note that body position in A Lanzada included stretched arms; López- Costas, 2015;seealsoFig. 1). Interestingly, here we find thatL07.20 (radius over the thorax at the ~L1 vertebra) and L06.22 (humerus) had higher mercury concentrations than the average of bone type 2. Therefore, bone could be acting as a source of mercury to the soil, but also as a sink for the mercury released from the body's soft tissues during decomposition (cf. Álvarez-Fernández et al. (2021)). 4.2. Bone structure indices Bone IR indices related to hydroxyapatite mineral structure (IRSF and MMI) correlated negatively with mercury content, suggesting that mercury tends to be higher in newly-formed bone tissue, and lower in mature bone tissue. The relationship with the mineral maturity index (MMI) and the infrared splitting factor (IRSF) indicates that Hg content increased in bone areas that were recently created or remodelled (Farlay et al., 2010). Hydroxyapatite content and crystal size increase with bone maturity, negatively affecting Hg retention in the bone structure since large crystals have lower specific reactive surface for binding. As expected, carbonate/ phosphate (CP) and amide I/phosphate (AmP) ratios correlated positively with bone mercury. This is consistent with their negative correlation with IRSF and MMI. New-formed and remodelled bone tissues have a higher proportion of collagen fibres and higher rates of C0 3 /PO 4 (Akkus et al., 2003) since they are less mineralized. Differences observed between new and remodelled bone (high mercury) and mature bone (low mercury) could respond to two different processes: i) an ante-mortem process: mercury is released or diluted due to bone tissues maturation and mineralization, since both large crystals and lower collagen content limit mercury retention; ii) a post-mortem process: mercury released during body decomposition is easily accumulated in new and remodelled bone due to its larger collagen content and smaller hydroxyapatite crystal size. IRSF and MMI showed no differences regarding individuals, consistent with previous research showing independence of the degree of bone mineralization with respect to the age and sex of adults, and more mature individuals (Boivin and Meunier, 2002). Differences in the degree of mineralization seem to be related to the age of the tissue (Bala et al., 2013). For N. Álvarez-Fernández et al. Science of the Total Environment 851 (2022) 158015 8 bone tissue, the degree of mineralization depends on the turnover rate (Boivin, 2007). However, differences were not found for IRSF and MMI between bone types 1 and 2, that have different turnover rate (Hedges et al., 2007). There are no differences in AmP between L06 and L07, in agreement with studies from biopsies (Danielsen et al., 1994) and autopsies (Wang et al., 2002), in which differences in collagen structure –in nonpathological individuals –were only found for individuals older than ~60 years. These studies indicated that changes in collagen are ageindependent and that changes in the turnover rate may exert some control. Pathologies like osteopenia can modify the bone turnover rate altering collagen content (Bala et al., 2013;Paschalis et al., 1997), which may explain why AmP is higher in L06 and L07 compared to L01, as L01 was an elderly woman with osteopenia and L06 and L07 were a juvenile and a mature adult with no pathological markers. AmP is also higher in bone type 1 than in type 2, perhaps also related to their different turnover rates (Akkus et al., 2003), and the fact that long bones are largely mineralized with a thick cortex, especially those which are related to locomotion. Consideration also may be given to the micro-sampling strategy, and to the fact that no metaphyseal areas were sampled for L06, so small differences between samples can be related to bone variability and only consistent correlations are considered here. 4.3. Mercury intra-skeletal variability The three models (G0, G1 and G2) explained an elevated proportion of mercury variability as well as that of bone hydroxyapatite (>90 %), while for bone collagen the percentage of variance explained was high to moderate(~60 %). G1 and G2 differed on their predictive powerfor soil mercury; G1 explained 17 %, whileG2 explained60 %. In G1, bone mercury had a no significant effect on soil mercury with a total path effect of 0.18. For G2, bone mercury had a significant positive effect on soil mercury with a total path effect of 0.37. The PLS-POS groups seemed to be responding to differences in bone ossification (cortical thickness) and surface chemical alteration (abrasion), being both more significant factors in G1 than in G2. Model G2 is most like G0, indicating that the samples in this segment dominated the general model. The main difference between the POS models relies in the bone component's interaction with bone and soil mercury. The bone mineral phase had a significant impact on bone Hg in the three models, while bone collagen was only relevant in G1. The bone mineral phase components affected soil mercury in opposite ways in G1 and G2. Thus, POS models seem to have captured intra-skeletal variability related to the weight of the bone components on mercury accumulation. Hydroxyapatite (HAP) had a positive effect on bone mercury in G0 and G2, while it was negative inG1. Collagen was not relevant in G0and G2 but dominated bone mercury variability in G1. Therefore, the PLS-POS groups responded to differences in the bone fraction controlling mercury content. In G1, the control was exerted by the organic fraction of the bone (i.e., collagen), while the inorganic fraction (i.e., hydroxyapatite) controlled the mercury content in G2. Bone is composed on average of ~70 % HAP, ~20 % collagen (type I), and ~8 % of water per weight (Augat and Schorlemmer, 2006;Currey, 2008). As bones of type 1 are expected to have higher turnover rates, then the proportion of HAP/ collagen will be presumably lower than in group 2 (Akkus et al., 2003). Thus, small increases in collagen in type 1 bones could have a larger impact on mercury accumulation, emphasizing the contribution of the less abundant component. This difference in the bone fraction controlling bone mercury content reflects a different behaviour related to the antemortem status of the bone tissue structure. Mercury tends to bind to soft bases like S, P, and N; bonds with S-containing groups being especially stable (Schuster, 1991). HAP can stabilise Hg as (Hg) 3 (PO 4 ) 2 (Cervini- Silva et al., 2021). Bone collagen structure is characterised by a repeating amino acid motif (Gly-X-Y), where X and Y can be any amino acid (Hulmes, 2008). All amino acids contain N-groups, but Cysteine (-SH, Cys) and Methionine (-S-CH 3 , Met) also contain S-groups. Gauza- Włodarczyk et al. (2017) estimated that the amount of Cys and Met per 100 g of collagen in bovine bone were respectively 0.32 g and 1.07 g. Hence, there is the potential to form Hg\\S compounds, and they may have an ante-orpost-mortem origin. To our knowledge, the mechanisms of mercury incorporation intothe bone during life are unknown, but this pathway cannot be discarded. On the other hand, during the main phase of soft tissues breakdown in the post-mortem span, the reducing conditions needed to form Hg\\S compounds are met (Janaway et al., 2009;Schuster, 1991). The presence of compounds that bind mercury strongly explains why the pathway from bone mercuryto soil is not significant when collagen controls bone mercury content, i.e., the release of mercury to the soil may be very low until the collagen structure is strongly degraded. The relationship between hydroxyapatite crystallinity (HAPc) (i.e., degree of crystallinity) and bone Hg was both negative and positive in the three models suggesting a dual role of bone crystallinity. This can be related to the potential of the skeleton to be both primary/secondary source and sink (Álvarez-Fernández et al., 2021). The models suggest that HAPc affects both bone Hg release (primary/secondary source) and retention/accumulation (sink; we cannot discern retention of bone ante-mortem content from mercury incorporated as the result of the interaction with soft tissues during their decomposition). In archaeological bone, as hydroxyapatite crystallinity increases, bone structure is altered (Nielsen- Marsh and Hedges, 2000) and, thus, as found here, a negative relationship with collagen and bone hydroxyapatite structure is expected. In G1 HAPc has a direct negative effect on bone mercury, as bone mercury is controlled by collagen. But there is also a slightly significant positive indirect effect, related to HAPc as a proxy of bone alteration. Despite bone mercury being controlled by collagen, it is very likely that HAP also contains mercury and as HAP is altered bone mercury will be lost. In G2, HAPc has a direct positive effect on bone mercury, reflecting that HAP controls mercury content in this group. When considering the ante-mortem condition of these bone samples, the higher the HAPc the higher the proportion of HAP (Akkus et al., 2003). But, when HAPc increases as a result of post-mortem bone alteration, the any accumulated mercury in HAP can be released. Location (i.e., the sample position regarding the pelvic area) has a negative effect both on bone and soil mercury content. Samples close to the thoracic/abdominal area had higher mercury concentrations compared to those farther away. This is consistent with what has been already discussed in Section 4.1., and related to the location of the main target organs for mercury in the thoracic/abdominal area (Berlin et al., 2015;García et al., 2001). It is known that bones may act as sink of some elements during body decomposition and then release them to the burial environment (Hedges, 2002). The relationship between location and bone mercury content strongly suggests that bones act as a temporary sink for mercury. 4.4. Role of soil components Most soil primarysilicates lack the capacity to retainmercury (Schuster, 1991), consistent with their dilution effect on bone and soil mercury content. The higher the amount of primary silicates the fewer soil components with the potential to retain mercury. Clay had a slightly negative effect on bone mercury content. The silt+clay fraction is a reactive component of the soil and tends to be enriched in metals (Qin et al., 2014). Therefore, a higher silt+clay content leads to an increase in mercury concentrations (Álvarez-Fernández et al., 2021;Cooke et al., 2020; Schuster, 1991) that limits the amount of mercury that remains available to be stabilised in bone.However, the weightof the clay in the threemodels is low (−0.13 G0, 0.01 G1, −0.17 G2). Bone mercury shows a significant direct positive effect on soil mercury in G0 and G2. This supports the idea that skeletons are potential sources of mercury to the soil/sediments in the burial environment (Álvarez-Fernández et al., 2021). It also supports the contention that bone mercury concentrations are not diagenetic when the geology of the area does not contain mercury ores, based on the data from archaeological skeletons (Emslie et al., 2015;Kepa et al., 2012; Rasmussen et al., 2013a;Walser et al., 2019;Yamada et al., 1995). This is also likely the case for the burials studied here because soil/sediment mercury concentrations are on average 0.9 ± 0.7 ng g −1 in the dune sands layer where the skeletons were buried (Álvarez-Fernández et al., N. Álvarez-Fernández et al. Science of the Total Environment 851 (2022) 158015 9