Physicochemical properties (colour, grain size, pH)
Samples collected outside the burials were lighter and had higher hue (dune sands, L*: 71.7±3.2; C*: 12.7±1.8; h: 80.9±1.9) than samples collected inside the burials (Necrosol, L*: 61.8±4.4; C*: 14.3±0.9; h: 79.1±0.9) (SI_Fig. 1). However, both had a larger yellow than red component (yellow, b*: 13.7±1.2; red, a*: 2.5±0.5). Results from grain size analysis indicates a predominance of sand fractions, being ~85% in Necrosol and more than 95% in samples outside the burial. In addition, Necrosol has higher content of gravel (2.76±1.7%), fine sands (22.22±1.5%) and silt+clay (9.79±2.81%) compared to samples outside the burials (gravel:1.46±2.14%; FS: 19.62±3.81%; SC: 2.32±1.62%) (SI_Fig. 1). pH values reflect alkaline conditions, although samples from dune sands had larger values (9.3±0.1 in water, 8.8±0.2 in KCl) than those of the Necrosol (8.9±0.2 in water, 8.3±0.2 in KCl) (SI_Fig. 1).
Both soil horizons of the paleosol (Ab-C) had different characteristics. The colour of the buried epipedon (Ab) was dark brown; being a*, b* and chromaticity higher (a*: 3.3±0.8; b*: 15.7±1.7; C*: 16.0±1.9) than in the dune sands (C horizon) (a*: 1.7±0.4; b*: 11.2±1.3; C*:11.4±1.3). Whereas luminosity and hue were lower in the Ab (L: 63.5±1.9 and h: 78.0±1.4) compared to the C horizon (L: 75.8±4.0 and h: 81.4±1.12) (SI_Fig. 1). Grain size showed a predominance of sands in both horizons (more than 70%), although medium sands were more abundant in the C horizon, while coarse and fine sands were more abundant in the Ab. Gravel content was three times higher in Ab (19.6±6.3%) than in C (6.0±9.7). The difference is much larger for the silt+clay: the content in Ab (6.4±3.3%) is about ten times higher than in C (0.6±0.8%) (SI_Fig. 1). pH results indicate higher alkalinity in C (9.4±0.1) than in Ab (8.9±0.1), being pH very homogeneous within each soil horizon.
Elemental composition (LOI, C, N and XRF)
The chemical composition is represented in SI_Figure 1. Necrosol’s LOI (0.68 ±0.15%) and N content (0.03±0.01%) were higher than that of dune sands (LOI: 0.37±0.07%; N: 0.014±0.006%) (SI_Fig. 1), while C showed the opposite (Necrosol: 3.36 ±0.23%; outside burials: 3.96±0.51%). Regarding the other elements, some of them (i.e., S and Si) showed high concentrations in punctual samples, while other (i.e., Al, K, or Cr) presented more homogeneous concentrations among samples. Silicon, Ca, Rb, Sr, Zr and U concentrations were higher in samples outside the burials. In contrast, P, Cu, Zn and Br concentrations were more elevated in Necrosol samples.
In the paleosol, LOI values were higher in the Ab (2.68±0.57%) than in the C horizon (1.09±0.37%). Nitrogen showed a similar distribution (Ab: 0.06±0.02%; C: 0.022±0.015%) (SI_Fig. 1) but carbon was lower in the Ab (1.7±0.11%) than in the C horizon (4.19±0.66%). The Ab also presented higher concentrations for Fe, Ti, Ga, Rb, Y, Pb, Th and Br, while the C horizon had higher S, Ca and Sr content. Phosphorus content was higher and Mn content lower in the Necrosol than in the Ab of the paleosol.
Spectroscopic analysis (FTIR-ATR)
The average and standard deviation spectra, as well as the average spectrum of the second derivative are represented in Figure 1. Six main absorbance areas can be observed: 3700-3400 cm-1, 2520-2510 cm-1, 1560-1300 cm-1, 1220-620 cm-1, and <550 cm-1. The standard deviation spectra shows that variability between samples is largest in the 1560-1300 cm-1 region, very high in the regions 1050-850 cm-1 and 600-500 cm-1, and moderate to low in the regions 1200-1050 cm-1 and 3700-3600 cm-1. The second derivative spectrum (Figure 1) enables to identified characteristic absorbances of soil components: quartz (1165, 1094, 1080, 798, 777, 693, 460 cm-1), K-feldspar (647, 535, 417 cm-1), carbonates (both calcite and aragonite, 2514, 1478-1411, 874, 859, and 712 cm-1), and clay minerals (i.e., kaolinite, 3694, 3668, 3647 and 3621, 1030, 1005, 911 cm-1); small amounts of mica (1005, 960, 527 cm-1) are also probable42,43,44. Very low absorbances around 3000-2800 and 1700-1600 cm-1 may correspond with low amounts of soil organic matter (SOM)42,45,46. While most of the vibrations in the region 1200-1050 cm-1 correspond to absorbances of silicate minerals, the standard deviation spectrum enables to identify a shoulder of moderate variability at 1200-1100 cm-1 that can be associated to vibrations of biogenic silica47,48,49,50,51.
We computed difference spectra52 by subtracting the average spectrum of the dune sands samples to the average spectra of the Ab samples, the samples collected outside the burials and those from T1 and T5 Necrosol (Figure 2). Negative differences are observed in the regions 2520-2510 (peaking at 2514 cm-1), 1560-1300 (peaking at 1411 cm-1), 900-850 (peaking at 874 and 859 cm-1) and 720-710 cm-1 (peaking at 712 cm-1). Positive differences are found in the regions 3700-3400, 1400-900, 600-500 and 480-420 cm-1. Negative differences correspond to carbonates’ vibrations, while positive differences correspond to quartz, clay and other silicates – and possibly also to SOM (OH vibrations around 3400 cm-1). Samples outside the burials do not show positive differences in the clay region (Figure 2) and only low negative values for the carbonates’ region. The other samples showed an increasing trend in both negative values and positive values following the sequence: C horizon outside burials à inside T5 à inside T1 à Ab horizon (Figure 2).
Principal components analysis
We performed a PCA using all analytical data (colour parameters, grain size, soil reaction, elemental composition, and selected absorbances of the IR data corresponding to soil components) obtained for the samples. The loadings of the 74 individual variables are in SI_Table 1. Ten components were needed to account for 90% of the total variance of the dataset, but only the first 5 contained a significant proportion of the variance of more than one variable and are the ones described here.
The first component, Cp1, accounted for 48.2% of the variance and showed large positive (>0.7) loadings for kaolinite (3694, 3651, 3619, 1029, 1005, 911, 693, 647, 606, 585, 531, 423 cm-1) and OH (3424 and 3215 cm-1) vibrations, chromaticity and colour components (a* and b*), total SOM indicators (LOI, N), organically-bound elements (Br), and metal elements (Fe, Pb, Th, Ti) (SI_Table 1). Moderate (0.3-0.7) positive loadings were also shown by the silt+clay fraction, silicate IR absorbances (431 cm-1), SOM absorbances (aliphatic SOM: 2922, 2879, 2853, and 2842 cm-1), and some major (Al) and metallic elements (Rb, Y, Ga, Cu, and Zn) (SI_Table 1). Variables with large (<-0.7) negative loadings include carbonates’ absorbances (1478, 1448, 1411, 874, 859, 712 cm-1), biogenic silica absorbances (1247, 1204, 1165, 1142, 1114 cm-1), hue (h) and luminosity (L*), soil reaction (pHw, pHk), total C, Ca, S and Sr, and medium and fine sands (SI_Table 1). Moderate (-0.43 to -0.65) negative loadings were also found for quartz absorbances (1081 and 1094 cm-1) and U.
Cp2 accounted for 9.9 % of the total variance, coarse sand and Zr showed a large positive loading, while P had a large negative loading (SI_Table 1). Moderate positive loadings were also found for many elements (Y, Rb, Mn, Si, Fe, Nb, N, Br, Pb), carbonates (i.e. aragonite, 859 and 712 cm-1), L*, soil reaction (pHw), and aliphatic SOM (2922 cm-1). Moderate negative loadings were found for fine soil components (silt+clay), silicates (quartz and kaolinite, 532, 449, 431 and 423 cm-1), Ca and Cu (SI_Table 1).
Cp3 (9.3% of the total variance) is dominated by large to moderate positive loadings of quartz (1094, 1081, 798, 777, 462, 431 cm-1), aliphatic SOM (2922, 2879, 2853, 2842 cm-1), biogenic silica (1165, 1142, 1114 cm-1), K and Rb. Carbonates (i.e., calcite and aragonite, 1411, 874, 859, 712 cm-1), clays (i.e., kaolinite, 3694, 911, 647, 606, 585, 531 cm-1) and S, have moderate negative loadings (SI_Table 1).
Cp4 (3.5% of total variance) shows no large loadings. Moderate positive loadings were obtained for Al, Ti, K, Si, U and colour components (a*, b* and chromaticity) and moderate negative loadings for Cr, Zn and hue (SI_Table_1). Cp5 (2% of total variance) is dominated by the anti-covariation of Cr, Si and U against Zn and Fe (SI_Table 1).
Figure 3 shows the components’ scores of the samples of the burials (outside and inside), and the paleosol. Positive Cp1 scores are found for the Ab and most of the samples from inside T1. The C horizon and samples collected outside the burials show negative loadings. Regarding Cp2, all samples from the paleosol have positive scores while almost all samples collected inside and outside the burials show negative loadings, exception made of a few samples outside the burials (Figure 3). Cp3 scores show an increase in value from the C horizon to the buried epipedon, to decrease again at the top of this horizon; samples from Necrosol show no clear trend, although those collected in the skulls have large negative scores while the rest have positive or slightly negative scores. No pattern is observed for Cp4 and Cp5 scores, both in the paleosol (Ab and C horizon) and burial samples (Figure 3).