Multi-isotope evidence for the emergence of cultural alterity in Late Neolithic Europe

See allHide authors and affiliations

Science Advances  22 Jan 2020:
Vol. 6, no. 4, eaay2169
DOI: 10.1126/sciadv.aay2169


The coexistence of cultural identities and their interaction is a fundamental topic of social sciences that is not easily addressed in prehistory. Differences in mortuary treatment can help approach this issue. Here, we present a multi-isotope study to track both diet and mobility through the life histories of 32 broadly coeval Late Neolithic individuals interred in caves and in megalithic graves of a restricted region of northern Iberia. The results show significant differences in infant- and child-rearing practices, in subsistence strategies, and in landscape use between burial locations. From this, we posit that the presence of communities with distinct lifestyles and cultural backgrounds is a primary reason for Late Neolithic variability in burial location in Western Europe and provides evidence of an early “them and us” scenario. We argue that this differentiation could have played a role in the building of lasting structures of socioeconomic inequality and, occasionally, violent conflict.


Prehistoric populations are frequently assumed to be internally homogeneous and bounded entities defined by their cultural distinctiveness, especially at regional scales of analysis, recalling ideas of traditional culture-historical archaeology that are being increasingly critiqued (1). The contexts in which, and processes by which, differentiation emerged in the past thus provide an important counter to this view, although one presenting considerable challenges in terms of archaeological recognition [e.g., (2)]. Even when diversity is explicitly sought within an archaeological culture, the results rarely show more than an “inconsistent commonality” [e.g., (3)]. One area that has long been of interest in this regard is variation in mortuary practices, an example of which is the nature of the relationship between funerary caves and megalithic graves in Neolithic Europe. Early interpretations suggested that variability in burial location might reflect change over time and/or different regional mortuary practices. However, it is now becoming clear that the use of funerary caves and of megalithic graves overlap both chronologically and spatially in a number of regions (e.g., Britain, Ireland, France, and Iberia) (46). Some accounts have argued for the blurring of the categories of built and natural burial places, particularly proposing their use at different stages of a multiphase funerary program spread across the landscape (6, 7). However, primary burial seems to predominate in both caves and megalithic graves across Western Europe [e.g., (810)]. Thus, approaches have generally relied on purported differences in “invested energy” with those interred in megalithic tombs being drawn from a more privileged segment of the community compared to those in caves (1113). This hypothesis might be supported by the finding of meaningful demographic differences between funerary locations that might broadly reflect a greater representation of adult males in megalithic graves (4, 10) and of differences in nitrogen isotope (δ15N) composition with megalithic people showing higher mean values (4, 14, 15), which has been essentially understood as reflecting higher lifetime amounts of animal protein. The alternative explanation is that both groups actually belonged to culturally different communities performing distinct funerary practices. Biodistance studies on dental traits have recently provided unexpected evidence of differences between some megalithic groups and nearby contemporary communities interred in caves (16). While DNA analyses demonstrate haplogroup diversity and identify independent familial groups in caves and megalithic graves (17), they have not been undertaken at sufficient resolution with these particular questions in mind.

The Rioja Alavesa region (north-central Iberia) has recently provided evidence of small but significant differences in bone collagen carbon isotope values (δ13Cbcol) between individuals deposited in caves [800 to 900 m above sea level (masl)] and in megalithic graves (550 to 700 masl) (18). This finding is robust because of the large sample size (n = 155), the restricted geographical area in which the sites are located, of some 250 km2 with no sites being more than 4 to 6.5 km apart (Fig. 1), and the strict chronological control. In Rioja Alavesa, caves were exclusively used as burial sites during the Late Neolithic (3500 to 2900 cal. BCE), while the use of megalithic tombs extends both earlier and later by some centuries (19). Consequently, only megalithic individuals that were considered contemporary (by stratigraphic context—e.g., discrete burial layers—and/or by direct radiocarbon dating) to those deposited in caves were included in the study. Such a spatiotemporal setting (essential to avoid confounding environmental and diachronic variables) is exceptional in the European record, where few opportunities exist to make a detailed comparison between those buried in megalithic and nonmegalithic graves. Evidence of regular interpersonal conflict (20) makes the region even a better “laboratory” for the study of the entangled socioeconomic, cultural, and, perhaps, ethnic relations grounded in communities performing different burial practices. The differences identified in δ13Cbcol values appear to be correlated with environmental factors (i.e., elevation, which plays a role in temperature, precipitation, and forest density, all of which can affect δ13C values in vegetation) (21), suggesting that use of the landscape was being divided at a very local scale. The reasons for this partitioning may involve differential social status (e.g., those interred in caves may be of lower standing with more restricted access to the valley’s arable fields) and/or economic specialization (e.g., upland herding versus valley farming) within the same community or, alternatively, culturally and/or biologically different populations using different burial locations and following distinct subsistence strategies (18).

Fig. 1 Location of the Late Neolithic burial sites of the Rioja Alavesa region and summed probability distributions of the available Late Neolithic/Early Chalcolithic radiocarbon dates from the funerary sites under study [(19)], pooled by burial location.

Locations (top): 1, Las Yurdinas II; 2, Los Husos I; 3, Los Husos II; 4, Peña Larga; 5, La Peña de Maranón; 6, La Cascaja; 7, El Montecillo; 8, Layaza; 9, El Sotillo; 10, San Martín; 11, Alto de la Huesera; 12, Chabola de la Hechicera; 13, El Encinal; 14, Los Llanos; 15, Longar. The radiocarbon dates (bottom) are modeled using OxCal 4.2.2.

Our knowledge of prehistoric settlement and subsistence practices in the Rioja Alavesa is limited, because of loss through historic agricultural activities and development and the lack of systematic survey. Regional palynological data for the period show depletion of forest and the presence of arable fields and pastures, consistent with an agropastoral regime (22), with Late Neolithic sheep/goat pens in some upland rock shelters (23) confirming a pastoral element. However, archaeobotanical and zooarchaeological evidence in funerary sites is too scarce to assess to what extent those communities burying their dead in caves and megalithic graves relied on farming and herding. The only evidence for crop cultivation comes from the presence of more than 3% cereal pollen in some cores [e.g., (22)]. Animal husbandry is supported by the finding of a few bones of domestic ungulates mainly comprising sheep/goat, cattle, and pigs. Hunting and gathering may have also played a role in subsistence, as suggested by the occasional finding of bones of red deer and wild boar and some nut shells (18). Given these resources, stable isotope analysis on human bone collagen alone cannot address potential differences in the proportion of foods consumed among and between burial locations. Lipid analysis of pottery may provide insights, but so far, there are no studies for the region. With this evidence in hand, little is therefore known concerning socioeconomic organization and hierarchy (24), let alone their potential relationship to subsistence.

The aim of this paper is to choose between the aforementioned alternatives through a multi-isotope investigation of the life histories of individuals interred in caves and in megalithic graves. To this end, two questions are posed: (i) At what age does the isotopic difference previously observed in δ13C of bone collagen appear and (ii) does this finding reflect different dietary practices or a division in the use of the landscape for the communities interring their dead in caves and megaliths?

The age at which the divergence in δ13Cbcol values between caves and megalithic graves appears has crucial interpretative implications. If seen after a certain age (e.g., adolescence), then it could relate to an increased differentiation within the same community (perhaps as a result of an individual’s acquired position or economic specialization). Conversely, if it is apparent from infancy or early childhood, then it could be linked to the existence of different communities following distinct subsistence economies and funerary practices within a partitioned landscape. To test this hypothesis, we measured dentine collagen carbon (δ13Cdcol) and nitrogen (δ15Ndcol) from 632 sequential samples extracted from first permanent molars (M1s) and second permanent molars (M2s) to track juvenile dietary life histories in 27 individuals from caves (Las Yurdinas II, Los Husos I, and Peña Larga; n = 14) and megalithic graves (Alto de la Huesera, Chabola de la Hechicera, and Longar; n = 13) of the Rioja Alavesa region (Fig. 1). Here, δ15N measurements are primarily used to establish the trophic position of an individual in the food web (25) and to assess possible differences in breastfeeding and weaning patterns (26). While 27 individuals may seem a small sample size, it is important to emphasize that the difference in δ13C between caves and monuments has already been demonstrated with a large number of adult and adolescent individuals (18). We are seeking the reason behind this differentiation.

Subtle isotopic differences may emerge through consuming different foods and/or through consuming the same foods but from different parts of the landscape, as represented in the Rioja Alavesa by uplands versus the valley bottom. The analysis of contemporary fauna could provide more detail on potential differences in subsistence practices, but our ability to produce a robust faunal baseline for the region is limited by the paucity of recovered faunal remains for the period (18). Taking another approach, enamel apatite carbon isotope values (δ13Cap) from the M1 and M2 of the aforementioned 27 individuals and from five additional individuals (one from Las Yurdinas II, one from Peña Larga, and three from Longar) are measured to estimate the dietary importance of plants through collagen-apatite spacing [e.g., (27)]. Strontium (87Sr/86Sr) and oxygen (δ18Oc) isotope measurements and strontium ([Sr]) and calcium ([Ca]) concentrations of enamel from these 32 individuals are analyzed in M2s. The 87Sr/86Sr ratios measured in human enamel derive from a combination of all foods and water consumed, which, in turn, reflect a complex mixture of weathered sediments, stream waters, precipitation, and prehistoric anthropogenic inputs in agricultural soils (28), while δ18Oc values provide information about sources of drinking water (29). Both values reflect an individual’s childhood subsistence landscape rather than an exact geological location (30). As a result of the process of biopurification, [Sr] values not only can be used as a proxy for trophic level (i.e., they decrease significantly from plants to animals) but can also reflect the use of Sr-rich food additives such as salt and underlying geological variability (31). As Sr substitutes for Ca in skeletal mineral, [Ca] values can also provide useful insights into these issues. [Sr]/[Ca] ratios are primarily used here to scale Sr inputs (32). These bioaveraging effects provide an opportunity to explore prehistoric group differentiation. Modern vegetation, water, and salt were analyzed to characterize the isotope compositions representative of the different geological formations that constitute the regional “isoscape” (33).


Dentine collagen carbon (δ13Cdcol) and nitrogen (δ15Ndcol) isotope data

A total of 632 dentine samples obtained from the M1 and M2 of 27 individuals met established quality criteria for well-preserved collagen (see table S1). Statistical comparisons between burial locations using a mixed-model nested analysis of variance (ANOVA) reveal significant differences in δ13Cdcol values, with cave individuals showing lower values (0.3‰ on average) than megalithic grave individuals [F1,24.94 = 5.22, P = 0.031, Cohen’s d = 0.881, statistical power (1-β) = 0.95] (Fig. 2 and fig. S1). When examining the values by age categories (0 to 2.9, 3 to 4.9, 5 to 6.9, 7 to 9.9, 10 to 11.9, and 12 to 16) calculated by equating formation stage with tooth age development (34), differences are seen from age 5 onward (table S2), except for the 10- to 11.9-year age category. Analyzing the results by sex, the divergence between caves and megalithic graves is only significant in males for the 3- to 4.9-year age category (table S3). No differences are seen in δ13Cdcol values between sexes either at a general level or within each burial location. There are no significant differences for δ15Ndcol whatever the comparison, a result consistent with the previously reported bone collagen values (18).

Fig. 2 Comparison between dentine collagen carbon (δ13Cdcol) and nitrogen (δ15Ndcol) isotope profiles of M1s and M2s from caves (green) and monuments (blue).

Shaded areas represent the mean values ±1 pooled SDs (√GV) at different age periods (0 to 2.9, 3 to 4.9, 5 to 6.9, 7 to 9.9, 10 to 11.9, and 12 to 16). Mean values obtained on the bone collagen (δ13Cbcol and δ15Nbcol) of cave and megalithic individuals at adulthood (18) are shown for reference. A schematic graph has been included to illustrate the correspondence between tooth developmental ages and anatomy (34) used to assign age to each dentine sample (where DEJ is dentine enamel junction, CEJ is cementum enamel junction, and VPDB is Vienna Pee Dee belemnite).

Observations pertaining to exclusive breastfeeding practice are limited by the low number of samples attributable to the first months of life, as a result of dental wear on the cusps of some M1s (table S4). Nevertheless, statistically significant differences are seen in the estimated age of complete cessation of breast milk consumption between those from caves (3.3 ± 0.7 years) and from megalithic graves (3.9 ± 0.5 years) (t = 2.332, df = 24, P = 0.028; Fig. 3 and cf. table S4). Although microdentinal samples represent an averaged (although relatively short) time span rather than a specific point in time, the estimated mean age is used for the sake of simplicity. No differences are seen in this regard between males and females at each burial location.

Fig. 3 Dispersion of weaning cessation ages by site type and by sex and site type.

Left, by site type; right, by sex and site type.

Enamel apatite carbon (δ13Cap) and carbonate oxygen (δ18Oc) isotope data

Sixty-four δ13Cap and δ18Oc measurements were obtained from the M1 and M2 of 32 individuals (table S5 and fig. S2). No differences are observed between burial locations in bulk δ13Cap values of M1 crown enamel, which forms between ca. 0 and 3 years (34). By contrast, those of M2 crown enamel, formed between ca. 2.5 and 8 years, show statistically significant differences between caves (−13.9 ± 0.8‰) and megalithic graves (−13.0 ± 0.7‰) (t = 3.246, df = 30, P = 0.003; Fig. 4), consistent with the pattern seen in δ13Cdcol values (cf. table S2). A pattern not reflected in δ13Cdcol is the existence of divergent M1 to M2 shifts in δ13Cap between burial locations, with a mean decrease of 0.3 ± 0.9‰ in caves and a mean increase of 0.6 ± 0.8‰ in megalithic graves (t = 3.108, df = 30, P = 0.004; fig. S3). This seems to be the reason behind the differences seen in the Δ13Cdcol-ap spacings of the M2 between burial locations (t = 2.300, df = 25, P = 0.030), with caves showing a smaller mean offset (6.0 ± 0.7‰) than megalithic graves (6.5 ± 0.6‰) (table S6 and fig. S4). No differences are seen between the sexes whatever the comparison.

Fig. 4 Dispersion of δ13Cap values by site type and by sex and site type.

Left, by site type; right, by sex and site type.

In contrast, the δ18Oc data show no clear difference by burial location at a general level, whether considering the M1 or the M2 values or the offsets between these (cf. table S5), despite the significant negative correlation identified between regional modern surface water δ18Odw values and elevation [correlation coefficient (r2) = 0.462, P = 0.003; rho = 0.647, P = 0.005; table S7 and fig. S5]. However, the data show a statistically significant difference between M2 values of males interred in caves (−4.7 ± 0.5‰) and in megaliths (−4.1 ± 0.6‰) (U = 12.0, Z = 2.309, P = 0.021) and also between males (−4.7 ± 0.5‰) and females in caves (−4.1 ± 0.6‰) (U = 8.0, Z = 2.315, P = 0.021) (fig. S6). The difference is also seen in the higher mean offset in males buried in caves between the M1 and M2 values (−1.3‰ versus ca. −0.6‰ in the rest of the sample).

Enamel strontium isotope (86Sr/87Sr) ratios and strontium ([Sr]) and calcium ([Ca]) concentrations

86Sr/87Sr and [Sr] and [Ca] values from bulk M2 enamel of the 32 individuals were used in combination with a series of environmental measurements to explore whether childhood locations and diets differed for the two sets of sites. Comparison of 87Sr/86Sr values between burial locations reveals significant differences (U = 54.0, Z = 2.798, P = 0.005; cf. table S5), with caves showing higher mean values (0.7087 ± 0.0008) than megalithic graves (0.7081 ± 0.0005) (fig. S7). This is consistent with their respective geologies (fig. S8), mainly consisting of Miocene and Oligocene sandstone in the valley and of Mesozoic limestone in the foothills, as supported by modern plant 87Sr/86Sr ratios obtained to map biologically available strontium (BASr) in the study area (table S8 and Fig. 5). There is significantly less variation among megalithic males compared to the rest of the population [SD = ±0.0003 (megalithic males), ±0.0007 (megalithic females), ±0.0008 (cave males), and ±0.0008 (cave females)] (Levene’s tests, P < 0.05; cf. fig. S7).

Fig. 5 Prehistoric human and modern plant strontium isotope values (87Sr/86Sr) plotted against latitude (ETRS89, 30°N).

Both the mean values estimated for local geological formations and the Rioja Alavesa region’s latitudinally cross-sectioned topography are shown for reference.

If the 87Sr/86Sr results are analyzed on a case-by-case basis to define “localness” (35), it can be seen that individuals whose values fall within the BASr range expected for a territory of less than 2 hours walk from each site predominate in both burial locations. There are small numbers of individuals whose values fall outside of the 2-hour catchments in both caves (2 of 16) and megalithic graves (3 of 16), although they still fall within the BASr range expected for the entire study region (tables S9 and S10 and fig. S9).

Comparison between burial locations reveals statistically significant differences in [Sr]/[Ca] ratios (U = 33.0, Z = 3.580, P < 0.001; fig. S10), which are mainly driven by the higher mean [Sr] values in megalithic graves [309.8 ± 163.3 parts per million (ppm)] compared to caves (127.8 ± 85.1 ppm), since [Ca] values do not differ between burial locations (cf. table S5 and fig. S2). This is consistent with the significantly higher [Sr] values obtained for both valley plants and salt (cf. tables S8 and S11). Assessing the relationship between [Sr] and other dietary proxies shows a significant correlation with δ15Ndcol values in caves (r2 = 0.421, P = 0.012; rho = 0.754, P = 0.002) and with δ13Cap values in megalithic graves (r2 = 0.364, P = 0.017; rho = 0.408, P = 0.131) (fig. S11). No patterns are seen when comparing [Sr] and either δ13Cdcol or δ18Oc values.

Comparing 87Sr/86Sr and [Sr] values shows less variability among those interred in megalithic graves, which may be evidence of a more restricted childhood subsistence landscape (i.e., the Rioja Alavesa lowlands) for all individuals but one (an outlier for both 87Sr/86Sr and [Sr], based on Z scores; Fig. 6). By contrast, greater variability is seen among those buried in caves, although only a few individuals overlap with those from the megaliths.

Fig. 6 Correlation between human strontium isotope values (87Sr/86Sr) and strontium concentrations ([Sr]).


While the existence of significant differences in bulk δ13C collagen values between Late Neolithic individuals interred in caves and in megalithic graves of the Rioja Alavesa region (n = 155) was established previously (18), its explanation remained unclear. Here, the analysis of fine-grained isotope data for the early lives of 32 individuals has served not only to confirm and assess the extent of these differences but also to clarify their potential origin. As collagen meeting the acceptance criteria cannot be diagenetic (36) and the rates of trace element leaching into enamel of the samples analyzed show no differential diagenesis between burial locations (37), the detected isotope differences must reflect in vivo practices.

Dentine collagen carbon isotope values (δ13Cdcol) (cf. Fig. 2 and fig. S1) extend to early childhood the differences previously identified in juvenile and adult bone collagen (δ13Cbcol) (18). This suggests that the relationship between both burial locations is neither the result of an individual’s acquired position nor of economic specialization within a single community, making use of both the valley and the mountain zones. Rather, distinct communities using different burial locations appear to have practiced different subsistence strategies and landscape use within a highly restricted region. Significant cave-megalith differences in δ13Cdcol values are observed across all age categories except for those younger than 5 years and those between 10 and 12 years (cf. table S2). Physiological processes and/or social age-related dietary differences may be behind these exceptions to what is otherwise a very robust trend. In the former case, they may be linked to varying enrichment in 13C (and 15N) for some individuals due to variability in breastfeeding duration and/or maternal-infant metabolic response (26). In the latter case, they may be principally attributed to the presence of δ13C (and δ 15N) dips around age 10 in the sequential dentine isotopic profiles of at least eight individuals (fig. S12), possibly as a result of the demands of rapid skeletal growth (38).

The δ15Ndcol results provide evidence for earlier weaning on average for those buried in caves than in megalithic graves (cf. table S4 and Fig. 3), indicating potentially distinct infant-rearing practices. This finding may be related to the significantly higher infant and adult female mortality rates in caves than in megalithic graves of the region (10), as the timing of weaning has strong implications for population dynamics and infant and maternal health (39). The duration of breastfeeding is usually shaped by religious beliefs, understandings of nutrition and health, child and parental identities, or practical considerations associated with the mother’s daily activities and mobility or sedentism. However, uneven intensities of unrest and food insecurity between the two groups might also result in distinct parental investment practices (40) or differentially increased physiological stress at early stages of life (38).

Enamel apatite carbon isotope values (δ13Cap) show significant differences in the M2 between funerary locations (cf. Fig. 4), consistent with the δ13Cdcol data, and suggest differences in complementary foods consumed during weaning and in childhood diet by those interred in caves and in megalithic graves. This finding is reinforced by divergences in the shifts between M1 and M2 values, which decrease in caves and increase in megaliths (cf. fig. S3). The fact that these shifts are not seen in δ13Cdcol values may reflect the preferential routing of dietary protein to collagen, which could mask the importance of plant carbon to the diet; apatite, by contrast, reflects whole diet (41). Differences in M2 Δ13Cdcol-ap values (cf. fig. S4) suggest greater plant consumption among megalithic children. Consistent with this, deciduous dentition shows significantly higher caries prevalence in megalithic graves than in caves of the region, which has been interpreted as a result of a more sticky texture and/or higher sucrose composition for childhood foods, such as gruel (42).

The significant negative correlation identified between regional modern water δ18Odw values and elevation, despite its subtlety, is consistent with the altitudinal effect in precipitation (43). However, this pattern is not seen in enamel carbonate oxygen isotope values (δ18Oc), which do not reveal any clear overall difference by burial location (cf. table S2 and fig. S5). This is consistent with the fact that (if atmospheric O2 is held constant) the δ18O in body water is influenced by both drinking water and dietary oxygen content. Moreover, physiological factors and/or culinary practices can cause variation in ingested water and food oxygen isotope values (44). Nevertheless, the data show a statistically significant difference between M2 values for males in caves and those for both males in megalithic graves and females in caves (fig. S6), which may be attributed to a different mobility pattern for males interred in caves perhaps due to their economic activity (e.g., herding, gathering, hunting, and trade).

The ratios between strontium and calcium concentrations ([Sr]/[Ca]) in those interred in megaliths are higher on average (cf. fig. S10), mainly due to their elevated [Sr] signals. A greater consumption of Sr-rich plant foods and/or salt (32) among individuals interred in megalithic graves might explain this difference. While the contribution of valley plants to megalithic diets seems clear and is the most obvious candidate, the consumption of salt cannot be excluded, especially in light of the presence of a number of hypersaline lakes in the Rioja Alavesa valley and of Neolithic salt production sites nearby (45). The correlations between [Sr] and δ15Ndcol in caves and δ13Cap in megalithic graves (cf. fig. S11) may tentatively be understood as status related, in the sense of differential consumption of higher trophic level animals (e.g., pigs and/or nursing herbivores) among those interred in caves with greater access to potentially more valued valley resources (e.g., cereals or salt, with the latter crucial for preserving meat) and of unequal access to different valley plots (e.g., open fields), plant products (e.g., bread and beer), or salt within those interred in megalithic graves.

The lower variability in strontium isotope values (87Sr/86Sr) among megalithic males compared to the rest of the population (cf. fig. S7) is consistent with a patrilocal residence system. It may even suggest inherited access to valley fields, as proposed for equivalent isotope findings in Early Neolithic Central Europe (28). Because megalithic graves in the Rioja Alavesa show a very regular pattern of distribution across the landscape, which has been attributed to a mixed economy focused on farming (46), substantial amounts (> 3%) of cereal pollen detected in megalithic graves (e.g., Longar) suggest the proximity of crop fields (47), and plant contribution to diet seems to have been significantly higher among those buried in megalithic tombs (this study), we propose that subsistence practices in megalithic populations principally relied on agriculture and, perhaps, seasonal exploitation of salt lakes, although, of course, some herding would also feature. By contrast, the comparison between 87Sr/86Sr and [Sr] (cf. Fig. 6) suggests a more diverse subsistence landscape and more restricted access to [Sr]-enriched valley resources among the communities using caves as burial places.

The differences in 87Sr/86Sr values by burial location are consistent with the different geologies in which the sites are located (cf. Fig. 5). This suggests that, despite some interaction between the groups (e.g., some 87Sr/86Sr results from caves are consistent with time spent in the valley, or foods consumed from there), a considerable degree of physical separation nevertheless existed between the two. Mobility and individual freedom can be curtailed by territorial tensions. Evidence for restricted landscape use and fixed territoriality may reflect both a preoccupation with protecting economic resources in a context of demographic pressure and territorial circumscription and reduced intergroup interaction due to perceived danger and insecurity (48). While the available data do not support the arrival of new people to the region, the chronological coincidence between evidence of population growth, the advent of the use of caves as burial places, and evidence of recurrent violence (20) suggests a pattern of conflict between culturally distinct groups.

All the observed differences in juvenile dietary patterns, subsistence practices, landscape use, and, possibly, residence systems between people buried in caves and in megalithic graves in the Rioja Alavesa region provide robust evidence for permanent socioeconomic and cultural (if not ethnic) differentiation linked to funerary practices and for lifestyle distinction at a fine-grained scale. These analyses permitted robust statistical tests, despite the sample size limitations and sites being no more than 4 to 6.5 km apart, which highlights the importance of pursuing new methodological avenues and scales of analysis to address long-standing archaeological questions. Comparable distinctions have been described in the Portuguese Estremadura based on dental affinities, implying a degree of actual biological distance (16). While not yet been supported by genetic studies at the local scale (17), we argue here that Late Neolithic funerary variability in the Rioja Alavesa, together with the observed isotopic trends, captures the active creation of cultural alterity in the Late Neolithic; the Early and Middle Neolithic periods, by contrast, show no comparable divisions in mortuary practices. Horizontal differentiation often leads to vertical differentiation (49), and our results suggest that this was the case in the Late Neolithic Rioja Alavesa. We further suggest that divergent beliefs and lifestyles could have been a source of tension and competition, on occasion escalating into lethal violence (20).


Sampling strategy

Lower M1s and M2s confidently assigned to Late Neolithic individuals by stratigraphic context and/or by direct radiocarbon dating were selected for analysis. Strict chronological control was the principal sampling criterion, particularly when confronted with multiple cumulative burials. Adult mandibles preserving both M1 and M2 and, at least, four more teeth were primarily targeted to guarantee replicability of the analysis and the preservation of material to permit future investigations. This criterion was established in agreement with museum curators and was a conditio sine qua non for sampling permission, as human teeth represent an important archive for a range of analyses. Molars that showed caries, severe dental wear, or damage (e.g., breakage and signs of burning) were excluded from sampling. Since their prevalence is generally low in the collections studied, the exclusion of these teeth is not thought to create any relevant biases in the final sample. From a potential sample of 119 individualized Late Neolithic adult mandibles available for the study region (18), only a subsample of 37 met all the criteria described above. Of them, five were excluded from the study for funding restrictions, seeking for parity in sample sizes between burial locations (caves versus megalithic graves) and between sexes. Sex estimation was performed following standard osteological criteria, based on craniomandibular and, where possible (i.e., in articulated skeletons), coxal traits (10, 18), corroborated by ancient DNA in three cases (LYII30, LYII41, and CH96) (17). While 32 individuals may seem a small sample size, the number of individuals that can be studied through high-resolution multi-isotope analyses is always limited because, despite continued methodological improvements, some techniques are still highly destructive. This is justified when performed to test explicit hypotheses, since these analyses provide nuanced dietary and mobility information for the early life histories of individuals that is not otherwise available. A distinct advantage of analyzing the M1 and M2 of adults is that they provide information about those who, by definition, survived into adulthood. The life histories provided by infant and child remains, by contrast, may be biased by the simple fact that they did not survive.

Modern plant, fresh water, salt water, and salt samples were also collected to characterize the regional isoscape. Branches and leaves from trees were sampled at different geological formations of the region, targeting diverse sampling locations within each formation (cf. table S8). Where possible, trees growing close to farmlands were avoided to prevent from potential contamination sources (e.g., pesticides) (35). Two samples were taken from each tree to check the validity of the results. Freshwater samples were collected from different locations along the course of several streams and rivers of the region. Saltwater samples were collected from the main local hypersaline lakes and from saltwater springs of nearby areas, where solid salts produced that there were also collected for analysis.

Dentine collagen carbon and nitrogen isotope analysis

The M1 and M2 of 27 individuals (14 from caves and 13 from megalithic graves) were analyzed. Teeth were recorded using high-resolution photogrammetry, shot-blasted with aluminum oxide pellets to remove surface debris, then embedded in Herculite II (a high-strength gypsum molding material), and sectioned longitudinally in half with a Buehler IsoMet low-speed diamond saw with a micrometer gauge, an abrasive wafering blade, and a cooling water bath (26). Collagen extraction was carried out following a modified University of Cape Town method (50) with sequential microsampling undertaken following a new method developed at the Research Laboratory for Archaeology and the History of Art (RLAHA), University of Oxford (40). One complete tooth half was demineralized in 0.5 M HCl solution at 4°C for 2 weeks and rinsed three times with deionized water until pH neutral. The demineralized tooth half was then cut longitudinally, and the posterior side was removed to minimize incremental growth overlapping as the pattern of dentine deposition follows an inward and rootward direction. The thin cementum layer was removed by mechanical cleaning of the tooth root surface, and any identified secondary and tertiary dentine was avoided by reaming out the pulp cavity before sampling. The remaining ~2-mm-wide longitudinal tooth slice, still retaining its original shape, was sampled sequentially from crown to apical root (mesial side) using a 1-mm-diameter KAI Medical biopsy punch with plunger, obtaining ca. 15 samples on average per tooth. This method is more anatomically sensitive than linear slicing. Resulting dentine microsamples were labeled according to a numerical sequence of dentine segments and their approximate age assigned by comparing the dentine sections’ anatomical location to the stage of dental development, following AlQahtani et al. (34). Demineralized dentine samples were freeze-dried before being weighed into tin capsules and loaded into a SERCON 20-22 continuous flow ratio mass spectrometer coupled with an elemental analyzer at the RLAHA.

Analytical precision was ±0.2‰ (1σ) for δ13Cdcol and δ15Ndcol based on repeated analyses of internal (alanine, marine seal, and cow) and international standards (International Atomic Energy Agency 600, caffeine). Dentine samples were generally not measured in duplicate, due to the small yields resulting from the sampling technique and to the impossibility of precisely replicating a sample elsewhere on the same tooth due to the complexity of dentine growth (51). However, a comparison was conducted on 28 samples with sufficient yields for duplicate measurement, resulting in a mean difference of 0.1 ± 0.1‰ (δ13Cdcol) and 0.2 ± 0.2‰ (δ15Ndcol) between the two runs. This is within measurement error and, moreover, the comparison between runs shows no significant differences (paired t test; δ13Cdcol: t = 0.539, df = 27, P = 0.595; δ15Ndcol: t = 0.211, df = 27, P = 0.834). Measurements of samples from caves and monuments were randomized to further ensure that no minor differences in the individual runs on the mass spectrometer would affect the results. Collagen preservation quality was checked according to standard criteria, including collagen yield, atomic weight C:N ratio, and %C and %N (52). With regard to atomic weight C:N ratios, in a conservative approach given the lack of duplicate measurements, a constrained C:N range of 2.9 to 3.4 was used. Last, only samples weighing ≥0.4 mg were considered reliable, since lower weights may affect mass spectrometer results [e.g., (53)].

Carbonate carbon, oxygen, and strontium isotope analyses and strontium and calcium concentrations

The M1 and M2 of 32 individuals (16 from caves and 16 from megalithic graves) were analyzed. M1s were only sampled for carbon and oxygen isotope analyses, while M2s were sampled for carbon, oxygen, and strontium isotope analyses, as well as strontium and calcium concentration measurements. The teeth were cleaned by abrasion with a diamond drill to remove the outer layer that could be contaminated with exogenous strontium. Once the outer layer was removed, further drilling was carried out to obtain 10 mg of powder for M1s and 50 mg for M2s. Care was taken to remain away from the dentine-enamel junction that is also more prone to diagenesis. If tooth fragments broke during drilling, then the surface of these fragments was removed by drilling. The remaining was then crushed with a mortar and pestle with the powder obtained by drilling. About 10 mg of the homogenized tooth enamel powder of the M2s was separated from the rest for carbon and oxygen isotope analyses.

The tooth enamel powder of both M1s and M2s was pretreated with 1 M buffered calcium acetate acetic acid (CH3COOH) solution in excess for 30 min (54). After treatment, samples were rinsed three times with MilliQ water and dried overnight at 50°C. About 1 mg of pretreated tooth enamel was then measured on a Nu Perspective Isotope Ratio Mass Spectrometer with a NuCarb carbonate preparation device at the Vrije Universiteit Brussel (VUB, Brussels, Belgium). Internal MAR2 (carbonate), ENF (enamel), and CBA (calcined bone standards), as well as Iso-Analytical IA-R022, were used. Repeatability of MAR2 measurements (n = 30) yielded an SD of reproducibility of 0.17 and 0.08‰ for δ13Cap and δ18Oc, respectively.

The remaining powder of M2 tooth enamel was pretreated with 0.1 M acetic acid before strontium extraction. Water samples were first evaporated on a hot plate at 100°C. The pretreated enamel powder, water deposits, and salts were then placed in Teflon beakers, dissolved in 14 M sub-boiled HNO3 for 4 hours at 110°C, and then left to dry overnight.

The plant samples were dried naturally, crushed with a coffee grinder before being ashed at 650°C in a muffle furnace by step heating for up to 8 hours. The ashed material was also placed in Teflon beakers and dissolved using a combination of 14 M sub-boiled HNO3, 23 M sub-boiled HF, and 10 M ultrapure HCl at 110°C. Once completely dissolved, the solutions were dried overnight.

The dried samples were redissolved in 2.5 ml of 2 M sub-boiled HNO3. A total of 0.5 ml was extracted for strontium concentration measurements (see below). The strontium of the remaining 2 ml was then extracted and purified following the protocol described in Snoeck et al. (55) and measured on a Nu Plasma multicollector-inductively coupled plasma Mass Spectrometer (Nu015 from Nu Instruments, Wrexham, UK) at the Université Libre de Bruxelles. During the course of this study, repeated measurements of the NBS987 standard solution yielded 87Sr/86Sr = 0.710239 ± 0.000026 (2 SD, n = 14), which is, for our purposes, sufficiently consistent with the mean value of 0.710252 ± 0.000013 (2 SD, n = 88) obtained by thermal ionization mass spectrometry instrumentation (56). All the data were corrected for mass fractionation by internal normalization to 86Sr/88Sr = 0.1194. All the sample measurements were normalized externally using a standard bracketing method with the recommended value of 87Sr/86Sr = 0.710248 for NBS987 (56). The 83Kr and 85Rb beam intensities are continuously monitored in routine to prevent any spectral interference; the intensities are always below 0.05 mV, which is negligible compared to the total Sr beam intensity (7 to 8 V). Procedural blanks were considered negligible [total Sr (V) of max 0.02 versus 7 to 8 V for analyses; i.e., ≈0.3%]. For each sample, the 87Sr/86Sr value was reported with a ±2σ error (absolute error value of the individual sample analysis—internal error).

The 0.5 ml of dissolved sample (see above) was used to measure strontium and calcium concentrations. Following dilution with 2% HNO3, Sr and Ca concentrations in the sample digests were determined using a Thermo Scientific Element 2 sector-field ICP mass spectrometer at the VUB (Belgium) in low (86Sr and 88Sr) and medium (43Ca and 44Ca) resolution using indium (In) as an internal standard and external calibration versus various reference materials (SRM1400, CCB01). On the basis of repeated digestion and measurement of these reference materials, the analytical precision (±1 SD) of the procedure outlined above was estimated to be better than 3%.

Modern water oxygen isotope analysis

Water samples were filtered after being collected and deposited into dark glass microtubes to avoid light exposure. They were analyzed using cavity ring-down water isotope spectroscopy in a Picarro L2130-i isotopic H2O analyzer at the RLAHA, University of Oxford. Each sample was generally analyzed three times (three consecutive replicate injections) alongside a set of two liquid water laboratory reference materials (LNW ZNW) that had previously been calibrated to the Vienna Standard Mean Ocean Water (VSMOW)-Standard Light Antarctic Precipitation international isotope scale. Analytical precision was ±0.3‰. All δ18Odw values were expressed relative to VSMOW (57).

Statistical tests

Shapiro-Wilk tests were used to assess whether the data were normally distributed. Levene’s tests were used to assess differences in the homogeneity of variances. For two-sample comparisons, Student’s t tests were used when the data did not depart significantly from a normal distribution; when they did, nonparametric Mann-Whitney U tests were used. Two-level mixed-model nested analysis of variance (ANOVA) tests were used to assess variation in δ13Cdcol value means among groups (i.e., caves versus megalithic graves) since the comparison involves one measurement variable (δ13Cdcol) and two nominal variables (individuals and burial location), the latter being nested (i.e., forming subgroups within groups) (58). Because of unequal group sample sizes, the Satterthwaite formula was applied (59). This allows the use of modified mean squares at each level, providing a better estimate of the effective degrees of freedom and, therefore, more accurate P values. Effect size was assessed using Cohen’s d (60). This statistic refers to the difference between two sample means expressed in pooled SDs, assuming approximately equal variances. Conventionally, a d value of less than 0.2 is seen as trivial, regardless of the outcome of the significance tests, while a value of 0.8 or more is termed large. 1-β was also calculated to measure the probability that a test is correctly rejecting the null hypothesis, especially when comparing small sample sizes. Pearson’s r/r2 and Spearman’s rho coefficients were used to assess correlations as appropriate. Z scores were calculated to detect the presence of outliers. A significance level of α = 0.05 and two-side testing was used for all the statistical tests. Means and SDs are used as measures for center and dispersion throughout the text, unless otherwise stated.


Supplementary material for this article is available at

Fig. S1. Density distribution of dentine collagen carbon (δ13Cdcol) isotope values, by site type and by sex and site type.

Fig. S2. Summary graph of dentine collagen nitrogen (δ15Ndcol) and carbon (δ13Cdcol) isotope values, enamel carbonate oxygen (δ18Oc) isotope values, enamel apatite carbon (δ13Cap) and strontium (87Sr/86Sr) isotope values, and strontium ([Sr]) and calcium ([Ca]) concentrations for the 32 individuals analyzed in the study.

Fig. S3. Dispersion of (δ13Cdcol) and (δ13Cap) ΔM1-M2 values for the 27 individuals selected for a life history approach.

Fig. S4. Dispersion of M1 and M2 Δ13Cdcol-ap values for the 27 individuals selected for a life history approach.

Fig. S5. Correlation between latitude and both M1 and M2 enamel carbonate (δ18Oc) and modern water (δ18Odw) oxygen isotope values.

Fig. S6. Dispersion of enamel carbonate oxygen isotope values (δ18Oc) by site type and by sex and site type.

Fig. S7. Dispersion of strontium isotope values (87Sr/86Sr) by site type and by sex and site type.

Fig. S8. Map with the geological formations of the Rioja Alavesa region and surrounding areas, with reference to the locations where modern plant, salt, and water samples were obtained.

Fig. S9. Strontium isotope ratios (87Sr/86Sr) for the individuals analyzed compared to the average BASr values for their burial sites and 15-, 30-, 60-, and 120-min walk distance catchments around these.

Fig. S10. Dispersion of [Sr](ppm)/[Ca](%) ratios by site type and by sex and site type.

Fig. S11. Correlation between M2s’ [Sr] values and (i) dentine collagen carbon isotope values (δ13Cdcol), (ii) dentine collagen nitrogen isotope values (δ15Ndcol), (iii) enamel apatite carbon isotope values (δ13Cap), and (iv) enamel carbonate oxygen isotope values (δ18Oc).

Fig. S12. δ13Cdcol (open purple squares) and δ15Ndcol (open orange triangles) isotope profiles of the M1 and M2 of the eight individuals showing clear dips around ages 9 to 11.

Table S1. Sequential dentine collagen carbon (δ13Cdcol) and nitrogen (δ15Ndcol) isotope values from the 27 individuals selected for the early life history approach.

Table S2. Comparison between cave and megalithic grave dentine collagen carbon (δ13Cdcol) and nitrogen (δ15Ndcol) isotope values.

Table S3. Comparison between cave and megalithic male and female dentine collagen carbon isotope values (δ13Cdcol).

Table S4. Estimated duration of exclusive breastfeeding and weaning process, age at complete weaning, and preweaning and postweaning dentine collagen nitrogen isotope values (δ15Ndcol) for the 27 individuals selected for the early life history approach.

Table S5. Summary table of adult bone collagen carbon (δ13Cbcol) and nitrogen (δ15Nbcol) isotope values, dentine collagen carbon (δ13Cdcol) and nitrogen (δ15Ndcol) isotope values (mean), enamel carbonate oxygen isotope values (δ18Oc), and enamel apatite carbon (δ13Cap) and strontium (87Sr/86Sr) isotope, and strontium ([Sr]) and calcium concentration ([Ca]) values for the 32 individuals analyzed in the study.

Table S6. Comparison between M1 and M2 crowns’ dentine collagen (δ13Ccol) and enamel apatite (δ13Cap) carbon isotope values, including Δ calculation, by site and burial location.

Table S7. Modern stream and lake water oxygen isotope values (δ18Odw) from the Rioja Alavesa region.

Table S8. Modern plant samples analyzed for the creation of the BASr map of the Rioja Alavesa region and strontium isotope (87Sr/86Sr) and strontium ([Sr]) and calcium concentration ([Ca]) values obtained from them.

Table S9. BASr (87Sr/86Sr ± 1 SD) for the local area (“local BASr”) and the average BASr values calculated for 15-, 30-, 60-, and 120-min walking distance catchments.

Table S10. Correspondence between the enamel apatite strontium isotope values (87Sr/86Sr) of the 32 individuals analyzed and the average BASr values calculated for the 0- to 60-min, 60- to 120-min, and >120-min walking areas from their respective burial locations.

Table S11. Strontium isotope (87Sr/86Sr) and strontium ([Sr]) and calcium concentration ([Ca]) values of salt and both spring and surface salt waters from the Rioja Alavesa region.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We are grateful to J. Fernández-Eraso, J. A. Mujika, J. Armendáriz, and N. Rivera for providing unpublished archaeological information; P. Ditchfield for mass spectrometry technical support; P. Claeys for granting access to the carbon and oxygen isotope laboratory; A. Weber for funding the isotopic H2O analyzer; J. Aguirre and J. Sesma for the facilities rendered for sampling human remains; and A. Plata and J. A. Abecia for allowing the collection of salt samples. We are grateful to the anonymous reviewers, whose comments have brought about a marked improvement in this paper. Funding: This research was funded by the John Fell OUP Research Fund [project title: Social differentiation and life histories in Late Neolithic northern Spain (EBD10940 - 151/102), to R.J.S., T.F.-C., and J.A.L.-T.] and supported by the Basque Government (POS_2015_2_0001, to T.F.-C.), the British Academy (Newton International Fellowship, NF170854, to T.F.-C.), the Research Foundation Flanders (to C.S.), and the Flemish Institute for Research and Innovation (IWT, to N.J.d.W.). Author contributions: T.F.-C., J.A.L.-T., and R.J.S. conceived the overall research. T.F.-C. and J.O. designed the sampling protocol and collected the samples. T.F.-C., C.S., N.J.d.W., A.C., and N.M. performed the laboratory work and generated the isotopic data. J.O. designed the Geographic Information System approach to characterizing BASr and to defining localness. T.F.-C. and J.O. conducted the statistical analyses and produced the graphic contents. T.F.-C., C.S., J.A.L.-T., and R.J.S. analyzed the results in combination with other archaeological data. T.F.-C., C.S., and R.J.S. wrote the manuscript with input from all authors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or in the Supplementary Materials. Additional data related to this paper may be requested from the authors. The human remains analyzed in the study are curated at the Museo Bibat de Álava (Basque Country) and the Museo de Navarra (Navarre), both in Spain.

Stay Connected to Science Advances

Navigate This Article