Environmentally adjusted δ13C thresholds for accurate detection of C4 plant consumption in Europe

Estimating the proportion of C₃ versus C₄ plants in human and animal diet is a key part of bioarchaeological, palaeontological and ecological research. Scholars worldwide have been investigating the spread of millet across Eurasia because this highly nutritious and drought-resistant C₄ plant can address various questions about past societies^1,2. This includes complex social structures, the adoption of new subsistence strategies, the adaptation to challenging climatic and environmental settings, as well as mobility and health status^3,4,5,6,7. In ecology and palaeontology, identifying C₃ and C₄ diets provides insight into (past) habitats, niche partitioning, and animal behaviours^8,9,10,11,12.

Stable carbon (C) isotope analyses—paired with nitrogen (N) isotope analyses when investigating collagen—represent the most efficient and preferred method to identify C₄ plant consumption from bioarchaeological and palaeontological skeletal remains. Due to different photosynthetic pathways between C₃ and C₄ plants, the ratio of ¹³C/¹²C isotopes (expressed as δ¹³C in ‰) is significantly more negative in C₃ plants (−35 to −23‰) compared to C₄ plants (−14 to −10‰)^13,14,15. Plant δ¹³C is transferred to the consumer’s body tissues along the food chain, following well-described and quantified fractionation processes^16,17, which leads to enriched δ¹³C values in the tissues of C₄ plant consumers compared to C₃ plant consumers. In general, bone or dentine collagen δ¹³C values above −18‰ and enamel or bone apatite δ¹³C values above −10‰ indicate a mixed diet of C₃ and C₄ plants^13,18,19. In contrast, δ¹³C values above −12‰ and −4‰, respectively, represent a pure C₄ diet^13,18,19.

However, specific environmental and climatic settings influence the plant’s δ¹³C value²⁰. For instance, C₃ plant δ¹³C values are more depleted under oceanic or Mediterranean climates, forest soil, dense canopy, elevated humidity or increased CO₂ concentration^20,21. Conversely, continental climate, aridity, salinity (including sea-spray effect), elevated temperature or high altitude tend to enrich plant δ¹³C values, which are the basis of the terrestrial food chain^20,21. The geographical location of the investigated site is also determinant. On average, skeletal tissues can be 1 to 2‰ lower in high latitudes compared to low latitudes in Europe^21,22,23. Although C₄ plants react differently to environmental and climatic factors^14,24, their δ¹³C values can vary across latitude as well^25,26. This implies that the threshold value for C₄ plant consumption needs to be adapted depending on the geographical location and environmental settings to avoid misinterpretations of body tissue isotopic composition. This paper fills this research gap to avoid an over- or underestimation of C₄ plants dietary intakes of premodern animals and human communities across Europe.

The main research questions addressed here are: (i) In which regions of Europe is it required to apply an environmentally adjusted δ¹³C threshold value for identifying C₃ versus C₄ plant consumption? (ii) What is the magnitude of this adjustment? (iii) Can we identify specific biogeographic parameters related to this isotope variability in Europe? This study draws on over 4000²⁷ published δ¹³C values from charred archaeological C₃ and C₄ grains derived from Isotope-Ratio Mass Spectrometry (IRMS)^{28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101}. We present an innovative and, to the best of our knowledge, unprecedented ecozone-based model framework that integrates multivariate environmental data to facilitate the identification of C₃ versus C₄ plant diets in bioarchaeological, palaeontological and ecological research.

Ecozone cluster model

Using multicomponent environmental datasets from topographical and climatic variables, we applied k-means cluster analysis to determine zones of similar environmental conditions across Europe. The model provides 20 spatial clusters including one cluster with unclassified values where not all conditions were equally met (Figs. 1 and  S1, Table 1). Due to numeric quantization of k observations, we labelled the clusters based on the percentiles of the respective data ranges using expressional combinations of temperature, climatic moisture index (CMI), and topography. Some of them are geographically restricted to specific regions, such as the cold and humid ecozones 4 and 5 in North-Eastern Europe, the mild and very humid ecozone 8 at the Atlantic coast, or the hot and arid ecozones in the south of Europe, northern Africa and the Near East (e.g., clusters 10, 16 and 18). Other ecozones are spatially more scattered, for example European high mountain ranges (e.g., clusters 11 and 15). The results can be compared to the biome-based ecoregions from the literature¹⁰² despite the reduced number of modelled ecozones.

The sites from which the investigated grains originate cover 15 out of 20 modelled ecozones (Fig. 1, Supplementary data 1), representing most of the European geographical and ecological diversity. However, their distribution is biased by past and modern human activity, archaeological sampling and analytical strategies. Not all ecozones are equally represented in the sample and the ecozones 3, 12, 14, 15 and 18 were excluded from the analyses due to small sample sizes.

Isotope diversity among grain species

Despite the isotopic diversity between the main C₃ grain species included in this study (Fig. S2A, B, Supplementary data 2), the two dominant crops of this sample, i.e. the Triticum (wheat, n = 1923) and Hordeum (barley, n = 1843) species, are largely overlapping in the northern and southern parts of Europe (Fig. S2B, Supplementary data 2). Differences from the Central/Western European samples are mostly caused by the wide geographical area represented by this region. In the UK, however, the two crops show notably distinct δ¹³C values (Fig. S2B, Supplementary data 2), which possibly reflects species-specific differences (e.g., the different timing of the vegetation period, which is earlier for barley, while wheat is more impacted by the summer conditions)^28,103 or the environmental diversity of the fields used to grow the different crops²⁹. Because human and animal diet is never based on one single crop, the C₃ grain sample was kept as one entity for the rest of the analyses. In contrast, the C₃ and C₄ grains show distinct δ¹³C values (Fig. S3C, Supplementary data 2), as expected from their different photosynthetic pathways^13,14,15. They are thus considered separately in the rest of the study.

Temporal isotope variability

Among the entire C₃ grains dataset, there is hardly any evolution of δ¹³C over time (linear model [lm] R² = 0.017, p = 1.75e-17; Pearson’s r value = 0.13, p = <2.2e-16; Fig. S3, Supplementary data 2). When distinguishing between the geographical subsets UK (lm R² = 0.044, p = 0.000276; Pearson’s r value = −0.21, p = 0.000276), Southern (lm R² = 0.03, p = 3.43e-17; Pearson’s r value = 0.17, p = <2.2e-16), Central/Western (lm R² = 0.215, p = 1.80e-32; Pearson’s r value = 0.46, p = <2.2e-16) and Northern Europe (including Denmark: lm R² = 0.078, p = 4.13e-18; Pearson’s r value = 0.28, p = <2.2e-16; and excluding Denmark: lm R² = 0.059, p = 1.55e-09; Pearson’s r value = 0.34, p = <2.2e-16), the positive correlation between C₃ grain δ¹³C values and the grain mean date is particularly weak (Fig. S4A–E, Supplementary data 2). In particular, the slightly stronger relationship observed for Central/Western Europe (Fig. S4B) is biased by the youngest samples from Central France, which exhibit particularly enriched δ¹³C values (Fig. S4C). At the ecozone level, a weak to moderate and significant increase in C₃ grains δ¹³C values can be observed for ecozone 1 (lm R² = 0.359, p = 5.39e-15; Pearson’s r value = 0.60, p = 5.39e-15) and ecozone 17 only (lm R² = 0.223, p = 0.00157; Pearson’s r value = 0.47, p = 1.57e-03) (Figs. 1 and  S5, Supplementary data 2). The C₄ grains dataset has a small sample size and each geographical area is represented by short chronologies, which does not enable any proper diachronic analysis (Fig. S6A). The slight decrease in C₄ δ¹³C values over time might thus be only considered statistically significant in Southern Europe despite the chronological gap of nearly a thousand years between the oldest and youngest cluster (Fig. S6B, Supplementary data 2). This implies that the C₃ and the C₄ grains datasets were not subdivided into different chronological phases for the subsequent analyses.

Geographical isotope variability

Splitting the C₃ grain dataset into geographical subsets (UK, Northern, Southern, and Central/Western Europe) shows that the median δ¹³C value of C₃ grains from Northern Europe is approximately 1‰ lower than that from Southern Europe (Fig. 2A, Table 2). This confirms the previous observations made on different types of samples such as faunal remains^22,23 and modern plants²¹. Yet it has to be stressed that the standard deviation (1 SD) is quite large for both regions (±1.35 and ±1.05, respectively), implying some overlap. Despite the high latitude, C₃ grains from the UK exhibit among the highest δ¹³C values across time (Figs. 2A and  S4), which can be related to the oceanic climate^22,23. In Denmark, the low δ¹³C values of the oldest half of the sample (c. 3700–3000 BCE) shift to particularly enriched δ¹³C values for the most recent half of the sample (c. 1000 BCE–1000 CE) (Fig. S7, Supplementary data 2). This might reflect changes in agricultural practices and soil management following its decrease in quality starting from the Neolithic period^30,31.

Using the whole dataset, we observe a significant decrease in C₃ grain δ¹³C values with increasing latitude (Fig. 2B, Supplementary data 2). From the median values calculated for each latitudinal bin (Table 2, Fig. 2B), the C₃ grains δ¹³C values from sites above 50° latitude are on average 0.54 to 1.72‰ lower than those of grains from sites at latitudes below 50° (Fig. 2B, Table 2). This confirms the mean offset of around 1–2‰ between Southern and Northern Europe. In comparison, there is a mean variation of 0.46‰ among the median δ¹³C values of the latitudinal bins above 50° and approximately 0.33‰ among the median δ¹³C values of the bins below 50°. The difference between southern and northern sites is therefore substantial, yet related to an increasing degree of variability towards the north. When excluding the UK and/or Denmark from this dataset due to the overall elevated values in these regions despite their northern latitude, the decrease in δ¹³C values with increasing latitude is accordingly even stronger and more significant (Pearson’s r value: −0.25 for the whole sample, −0.27 excluding UK and −0.30 excluding UK and Denmark, with a p-value < 2.2e-16 in each case; see Supplementary data 2).

At the modern country level, the C₃ grains from Lithuania (median −25.18 ± 1.16‰, n = 153) are on average nearly 2.5‰ lower than those from Jordan (median −22.86 ± 0.74‰, n = 46) (Fig. 3A, Table 2) which exceeds the previously defined offset of 1–2‰ between Southern and Northern Europe^21,22,23. On the contrary, the Jordan sample is on average 0.80‰ more enriched than those from Greece (median −23.50 ± 0.82‰, n = 383) or Italy (median −23.50 ± 0.97‰, n = 497) despite their shared southern location. Consequently, the North–South-dichotomy is not enough to characterize the different isotopic composition of grains from diverse parts of Europe and does not account for micro-regional environmental diversity. Moreover, the northernmost countries show standard deviations (SD) of 1.29‰ on average (1.12 to 1.67‰ in total), which is sensibly more than in most of the southern (from 0.40 to 1.43‰; mean: 0.92‰) and of the central/western countries (from 0.39 to 1.51‰, mean: 0.90‰) (see the ANOVA test in Supplementary data 2).

Similarly, the C₄ grain δ¹³C values from Lithuania (median: −10.83 ± 0.48‰, n = 20) are on average nearly 1‰ lower than those from Greece (median: −10.19 ± 0.15‰, n = 12), France (median: −10.15‰, n = 16) or Poland (median: −10.19 ± 0.15‰, n = 9) (Fig. 3B, C; Table 2). The sample sizes from Spain (n = 3, median: −10.69 ± 0.15‰) and from the Czech Republic (n = 1) are too low to be considered representative. This trend confirms previous studies from China with depleted C₄ grain δ¹³C values recorded at higher latitudes^25,26. The pattern is further supported by the larger C₄ grain δ¹³C dataset resulting from Alpha Magnetic Spectrometer (AMS) in Europe¹⁰⁴, showing generally lower δ¹³C values in Northern compared to Southern or Central Europe (Fig. S8). However, AMS stable isotopic data lack precision due to differences in calibration compared to IRMS and provide particularly wide and unusual δ¹³C ranges for C₄ grains¹⁰⁵. Therefore, these data cannot be used to extend the IRMS dataset in this study. Differences in local plant genotypes are considered more likely triggers for the isotopic variability than climate and water availability²⁴. Both C₃ and C₄ grains exhibit lower δ¹³C values in some northern regions relative to the rest of Europe, leading to regional variation in the δ¹³C threshold for identifying C₄ consumption.

Ecological isotopic variability

The geographical isotopic variability is related to environmental factors captured in the ecozone model. C₃ grain samples from Lithuania (n = 153), Estonia (n = 11), Finland (n = 44), and from parts of Denmark (n = 86) fall into the subhumid temperate lowlands of North-Eastern Europe represented by ecozone 5 (n = 294). Together with ecozone 17 (n = 42)—represented by grains from Norway only—these samples exhibit the lowest median δ¹³C values (−25.01 ± 1.25‰ and −25.08 ± 1.20‰, respectively) for charred C₃ grains in Europe (Fig. 4, Table 2). The high humidity, low temperature and low to moderately elevated topography of these ecozones can explain the depleted δ¹³C values²⁰. Ecozone 20 (n = 717), represented by balanced temperate plains scattered over Europe and including samples mostly from Denmark, northern Germany and England, exhibit a much higher median δ¹³C value (−23.00 ± 1.16‰) and its range hardly overlaps with the other northern samples. Beyond the influence of agricultural practices mentioned above^30,31, this can be explained by higher temperatures and moderate humidity characterizing this ecozone. In contrast, the 1007 samples from ecozone 17 and the 140 samples from ecozone 1 show the lowest median δ¹³C values among the sites below 50° latitude (−23.47 ± 0.82‰ and −23.35 ± 1.41‰, respectively). This reflects the mild and moist conditions of these transition zones at a mid-range altitude. In Southern Europe, the ecozones 7 and 13, representing warm highlands in the Mediterranean area, show the highest median δ¹³C values (−22.93 ± 1.35‰ and −22.56 ± 0.87‰, respectively), deriving from the drier and warmer climatic conditions. The ecozones 2, 8 and 19 are scattered over wide areas of Europe and are not related to extreme temperatures. Their isotopic ratios show intermediate values (Table 2). Ecozones 3, 12, 14, 15 and 18 cannot be included in this isotope investigation due to their small sample size.

Building on the substantial geographical and ecological variation in isotope values within C₃ (and to a lesser extent C₄) plants, it is essential to revise the commonly used δ¹³C threshold for identifying C₄ consumption, such as −18.0‰ for mammal collagen (for example, ref. ⁷). At each investigated site, the C₃ and C₄ grain δ¹³C values from this dataset (Fig. 5A) were used to create theoretical collagen δ¹³C values for a diet based exclusively on these crops (Fig. 5B)—which is no realistic diet for humans or animals and was only used for a first theoretical model. This resulted in site-specific estimations for an overall C₃ grain-based diet with 10% to 20% C₄ grain inputs (Fig. 5B and Supplementary data 1). Our model shows that at several sites from the Baltic and Nordic countries, human or animal collagen δ¹³C values below −19.0‰ (with a mean SD of 0.59‰)—and up to −19.68 ± 0.94‰ at Bėlis lake, Lithuania (n = 10), for example—already reflect a low C₄ input within a primarily C₃-grain-based diet. In the Mediterranean, the same C₄ input would result in collagen δ¹³C values above −17.0 ± 0.77 or −16.0 ± 0.57‰ (Fig. 5B and Supplementary data 1).

It is generally accepted that at least 20% of dietary protein from an alternative source (such as C₄ compared to C₃ crops) is required to be detected in collagen¹⁰⁶. However, with a model based on theoretical grain-based diets, a 20% C₄ input produces excessively elevated δ¹³C estimates for mammal collagen (Supplementary data 1). This highlights the limits of a model based on non-realistic grain-based diets, as despite the important role of grains in human diet, both humans and animals have more varied diets in reality—the latter even consuming mostly other parts of the plant and only rarely grains. And with collagen δ¹³C values reflecting the protein component of the diet¹⁶, plant intake contributes less to collagen isotopic composition than animal-derived proteins, which may alter this threshold value in omnivorous diets. Moreover, in regions where the diet includes significant proportions of mushrooms¹⁰⁷, forest-derived foods influenced by the canopy effect^108,109, and/or freshwater fish^110,111,112, consumers are exposed to more depleted δ¹³C values. These foods have a much stronger influence on collagen δ¹³C values than any enrichment from C₄ plants. And since such conditions are highly plausible in the Baltics and in Scandinavia, the threshold δ¹³C value for C₄ consumption might be even lower than those suggested by the model with 10% C₄ input (Figs. 5A, B and  6). On the contrary, potential marine food consumption¹⁷ or sea-spray effects in coastal areas¹¹³ need to be considered to avoid any over-estimation of C₄ plant intake.

A model involving all other food resources would go beyond the scope of this paper. But to accurately estimate the actual local threshold values for C₄ consumption, it is essential to use δ¹³C and δ¹⁵N values from as many local and contemporaneous food resources as possible, including crops (for example³²). Wild plants would further represent a better baseline for herbivore’s diets. Yet these are much more seldom in archaeological remains and their isotopic ratios are hardly represented in bioarchaeological studies so far. A model based on modern plants²¹ or on tree δ¹³C values¹¹⁴ can thus further serve as comparison δ¹³C baseline for herbivore’s diets. The combination of various isotope systems^113,115,116 and/or the application of mixing models¹¹⁷ are further powerful approaches to disentangle the diverse dietary sources.

Considerable isotopic variability is observed within each site, region and ecozone (Figs. 2–4 and 5A, Supplementary data 1–2), indicating that a range rather than a sharp threshold value is more appropriate (Figs. 5B and 6B). This isotopic variability is not only related to the gradual and complex variation of environmental settings across Europe, but is also intrinsic to the grains, as grain δ¹³C may vary for up to 0.5‰ within one ear despite same species and same growing conditions³³. Differing growing conditions, in particular various watering practices, can further impact local grain δ¹³C values by up to 1.7‰^24,28. Wild plants δ¹³C values may be good comparison references to disentangle anthropogenic and natural differences in water regimes^21,118,119. Another variability might be induced by diverging analytical uncertainty between datasets, as the data result from different laboratories and protocols and were obtained with different calibrations—which are relevant information to make sure that the datasets are comparable¹²⁰. Unfortunately, information on calibration, precision and accuracy of the published isotope data were mostly missing, thus it was not possible to assess the impact on the presented results. Yet when available, the error value for replicated samples was very low and still within the range of analytical errors.

An evolution of the grain δ¹³C values might be expected over time as well, especially when considering the various climatic phases that implied a great variation in moisture and temperature throughout Europe over the investigated time frame. The presented regional differences in grain δ¹³C values have therefore varied between the climatic phases (Fig. S10), yet the isotopic variations within each region over the various climatic phases is overall significantly weak (Supplementary data 2). A larger dataset for the most recent periods might reveal more effective trends. In this study, the impact of the chronological depth and imbalance covered by this dataset can be considered weak (Figs. S2–S5 and S10).

The isotopic difference reported in this study between regions and ecozone therefore remains strong enough to highlight environmentally-driven differences (Figs. 2–7, Table 2) even from grains that were possibly undergoing various cultivation practices. Our approach thus demonstrates that the threshold at −18.0‰ for consumer’s collagen is not universal and is mostly valid in the ecozone 2 (−17.51 ± 1.28‰), whereas the northern ecozone 5 and the widespread ecozone 17 require a threshold δ¹³C value closer to −19.0‰, i.e., −18.59 ± 1.12 and −18.67 ± 1.08‰, respectively. In high-temperature Mediterranean lowlands (ecozone 16), European plains (ecozone 19), and Atlantic regions (ecozone 20), the threshold shifts to −17.0‰ (i.e., −16.88 ± 0.90‰, −17.15 ± 0.74‰, and −16.71 ± 1.05‰, respectively; Figs. 6, 7, Tab. S1). In the arid southern regions (ecozone 13 [−16.31 ± 0.78‰], and partially ecozones 7 and 16), it approaches −16.0‰ (Fig. 7, Tab. S1). Yet in each case, the SD ranges from 0.74 to 1.26‰, stressing again the isotopic variability within ecozones.

To conclude, this paper offers both a European-wide δ¹³C baseline from archaeological charred C₃ grains and a threshold-value-model for environmentally adjusted identification of C₄ consumption from the site to the ecozone level across Europe (Fig. 7, Tab. S1). The grain baseline offers the advantage to consider a fundamental dietary resource (in particular for humans) as reference data and can be completed by local foods resources at the site level for more holistic and accurate interpretations. A comparison to wild proxies further enhances the results for animal diets. The threshold estimations are particularly suitable in bioarchaeological, ecological or palaeontological studies for which local plants or food resources are unavailable for calculating an isotopic baseline. However, it requires to account for the great degree of isotopic variability within each geographical entity as underlined by the SD values—a variability slightly increasing with latitude. In this context, the point-based approach (Fig. 5) provides more accurate yet geographically discrete data, while the interpolation-based approach (Fig. 6) offers ecologically sensitive estimates over large areas at the ecozone level, both related to some degree of uncertainty or variability. Future datasets could be used to test (and if necessary adjust) the interpolated values in areas of currently low site density. This innovative and context-sensitive ecozone clustering model based on temperature, humidity, and elevation thus enables more accurate interpretations of both animal ecologies and anthropogenic social and agricultural dynamics across Europe by avoiding over- or underestimation of C₄ consumption.

Isotopic dataset

The material used in this study consists of published δ¹³C values from charred grains derived from archaeological context, compiled into one single dataset²⁷. The data was collected from 75 publications until September 2025^{28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101}, also using the open access online repositories IsoArch¹²¹, MAIA¹²², Isotòpia¹²³, IsoMedIta¹²⁴ and CIMA¹²⁵. In total, this represents 4,210 δ¹³C values of C₃ and C₄ grains derived from 260 sites dated between 8000 BCE and 1800 CE. The represented C₃ plants are oat (Avena species, n = 58), rye (Secale species, n = 325), barley (Hordeum species, n = 1843), and wheat (Triticum species, n = 1923). Broomcorn millet (Panicum miliaceum, n = 57) and foxtail millet (Setaria italica, n = 4) represent the C₄ crops from this dataset. To facilitate visualization and pattern-recognition, the C₃ and C₄ grain datasets were considered separately for the various analyses due to their distinct δ¹³C values and due to the particularly small C₄ grain sample size. C₃ and C₄ grain δ¹³C values were then combined to address the question of identifying the introduction of C₄ crops in C₃ plants-based diets. Despite the fact that anthropogenic agricultural practices such as irrigation can impact grain δ¹³C values^24,28, altering the natural and ecological signal, crops remain an important proxy for human diet regardless of agricultural practices, as their isotopic composition would be transferred to human tissues all the same. The ecological differences are considered important enough across Europe to be detectable from grain isotopic composition despite anthropogenic alterations. Wild plants would have represented a better proxy for animal diet, however, this proxy is lacking from archaeological contexts—or could be derived from tree δ¹³C values for the most recent periods¹¹⁴.

Geographically, the research area spans modern Europe and the Mediterranean countries of the Near East, i.e., between 30 and 63° (N) latitude and between −8 and 45° (E) longitude. Yet the data is not evenly distributed and Denmark is over-represented in terms of sites, while Greece is over-represented in terms of number of samples. There are considerable gaps in several regions of Europe (see Fig. 1 and the isotopic dataset²⁷). Chronologically, this dataset covers most archaeological and historical periods and ranges from 8000 BCE to 1800 CE. In this context, it is important to stress that the oldest samples (8000–6000 BCE) exclusively originate from Greece, whereas samples from Northern Europe are predominantly younger than 1000 BCE—except for some larger site samples dated between 4000 BCE and 2000 BCE. Modern data created in the framework of experimental archaeology were not included in the dataset because of the controlled conditions in which they were produced and because of the different present-day atmospheric composition compared to pre-industrial periods¹²⁶.

Only data obtained from Isotope-Ratio Mass Spectrometry (IRMS) were selected for analyses. Despite the fact that these values are not comparable to IRMS isotopic values due to different calibrations¹⁰⁵, a European-wide dataset obtained from Accelerator Mass Spectrometry (AMS) in the context of radiocarbon measurements¹⁰⁴ was also used as a comparison dataset to verify whether the same trend is visible among both datasets. Importantly, grain δ¹³C values are sometimes (i.e., in 12 out of 64 publications used in this study) published in form of corrected values to consider the charring effect on the carbon isotope composition of archaeological C₃^127,128 and C₄ grains^129,130. But following the approach by Gron et al.³⁰, we are considering here only uncorrected values in order to enhance the comparability between datasets. This is considered to have no significant impact on this study’s results, as the charring effect is not systematic³³ particularly low on grain δ¹³C values (i.e., 0.06 to 0.18‰^{33,127,128,129,130} on C₃ and C₄ grains for a heat up to 300 °C) and remains below analytical errors for isotope analyses¹³¹.

Statistical analyses

All statistical analyses were performed using R software¹³² and the results and relevant values are summarised in Supplementary data 2. To determine the relationship between grain δ¹³C values and chronology or geographical location, we applied both Pearson’s correlation tests and linear models—the latter using the lme4 package¹³³. The Pearson correlation coefficient (Pearson’s r value) indicated the strength and direction of a linear relationship between two tested variables, the confidence interval gives the uncertainty range for the true correlation. While in the linear model, the R-squared values show the percentage of the dataset affected by the relationship and the p-values show the significance of the results. To determine the geographical scale at which significant changes in grain values occur, the dataset was divided into various bins. At the largest scale, the main regions of Europe, split into Northern Europe, Southern Europe, Central/Western Europe and the UK. This is not only convenient but also follows expectations based on previous research^21,22,23. At the smallest geographical scale, the dataset was binned according to the borders of modern countries. Boxplots were used for data visualization and one-way ANOVA tests were performed to test the difference in grain δ¹³C values between the investigated clusters. The results and relevant values for the ANOVA tests are summarised in Supplementary data 2). Because our dataset shows a particularly important isotopic variability, the results for each considered bin or entity/group are presented in the text using the median value (since the mean value is more sensitive to extreme outliers) and the related one standard deviation (1 SD). The tables are showing the mean, median, median absolute deviation (MAD) and 1 SD for each category/bin. In order to integrate an ecological dimension to the investigation of grain δ¹³C variability, the dataset was also binned into newly determined ecozones, as presented in the section below.

Environmental cluster

For the cluster analysis, we used these R-packages: terra¹³⁴, sf ^135,136, gtools¹³⁷, dplyr¹³⁸ and ggplot2¹³⁹, ggspatial¹⁴⁰, gridExtra¹⁴¹ for plotting, as described in the R-code provided in the open access repository to this paper¹⁴². We used an unsupervised k-means clustering approach based on three spatial predictors to differentiate environmental ecozones across Europe. Components include elevation (DEM), mean temperature (T), and the Climatic Moisture Index (CMI). A DEM derived from the USGS (United States Geological Survey, Global Multi-resolution Terrain Elevation Data 2010, https://earthexplorer.usgs.gov/; last accessed 19^th of June 2025). Monthly resolved climate variables for the period 1980-2018 were downloaded from CHELSA¹⁴³. The CMI represents a standardized water availability index, calculated as the ratio of precipitation (P) to potential evapotranspiration (PET) with CMI = P − PET. We used CHELSA v2.1 monthly CMI data based on the Penman-Monteith equation for PET and downscaled from ERA5 reanalysis. The grids were reprojected to a meter-based projection (EPSG:3857) and cropped to the extent of the geographic European landmass. The dataset (n = 468) were aggregated to a 1000 m resolution using bilinear interpolation prior to clustering. Monthly layers were averaged to create a multiannual mean (1980–2018). To allow comparison between variables with different scales and units, all raster layers were standardized using z-score normalization with mean and standard deviation (sd) of the raster (z = cell value − mean_raster / sd_raster). A regular grid of 1 km spacing was generated across the study area and the centroid of each grid cell was calculated for regular point sampling. To reduce edge effects and avoid NA values near the coast, centroids were restricted to land areas using a simplified buffer around the European landmass boundaries. At each centroid location, the values were extracted from the normalized rasters.

K-Means clustering

K-means clustering was performed using the extracted values, with the determined number of clusters k = 20. This dimensionality was chosen to balance regional ecological resolution with model interpretability, including NA values (replaced by −99999 during the cluster analysis to protect correct geographical raster reassignment). Each centroid was assigned a cluster label based on the combined environmental profile of elevation, temperature, and moisture availability. Cluster labels were spatially joined to centroid coordinates and rasterized back into a continuous spatial layer using the DEM grid as a template. The resulting map classifies observational sites into discrete clustered ecozones (Fig. 1).

To characterize each of the 20 ecozone clusters based on the data variability, we summarized the environmental properties of each cluster using the mean and standard deviation of the three input variables: Elevation (ELEV), T, and CMI. All values were normalized using z-scores prior to clustering, enabling direct comparison across variables. The clusters were then qualitatively interpreted based on their relative environmental signatures. For example, clusters with low temperatures, moderately dry conditions and high elevations were categorized as Cool|Moderately Dry|Alpine zones. Clusters with high temperatures and low moisture availability in low areas were described as Very Hot|Very Dry|Moderately Low. Intermediate clusters were labeled based on transitional or temperate climate conditions (Table 1).

Determining consumer’s δ
¹³C values

To determine the theoretically expected δ¹³C values of consumer tissues from a mixed C₃-C₄ diet, we first associated each C₃ grain δ¹³C value to a measured or assumed C₄ grain δ¹³C value. This means that for each site at which δ¹³C values of both C₃ and C₄ grains were available, each C₃ grain δ¹³C value got associated with the mean C₄ grain δ¹³C value of the site. However, most of the sites included in this study provided no C₄ grain. In this case, a theoretical C₄ grain δ¹³C value was determined for the site based on the observed values from this dataset. In most regions of Europe, the C₄ grain δ¹³C value is thus expected to be −10‰^32,72,87. Yet we observed that the C₄ grain δ¹³C values in Lithuania—and by extension presumably in the northernmost latitudes of Europe—were rather around −11‰⁶⁸. Similarly, C₄ grain δ¹³C values from the western Mediterranean area seem to be closer to −10.5‰^35,47. These regional values were thus used as theoretical C₄ grain δ¹³C values associated with the measured C₃ grain δ¹³C values at each site lacking C₄ grains (see summary in Supplementary data 1).

In a second step, 5‰ was added to each measured or theoretical grain δ¹³C value to mimic the fractionation offset that applies between the diet and the consumer’s collagenous tissues after consumption^16,17. We applied this to both C₃ and C₄ grains, resulting in theoretical end-members collagen δ¹³C values for 100% C₃ and 100% C₄ based diet, respectively. In a third step, we used these end-members values for each grain to create theoretical collagen δ¹³C values for a C₃-grain-based diet including either 10% or 20% C₄ grains. These three first steps were done at the grain level to minimize the loss of resolution and information when working with mean values. In a fourth step, we eventually calculated a mean δ¹³C value for these two types of diet at the site level (Supplementary data 1). It is fundamental to note that these fictitious diets based on 100% grains are not existing in nature and only represent a theoretical model using grains only.

Because a 100% grain-based diet does not exist in nature, the model using 10% of C₄ input was considered the most reliable basis for estimating the related collagen δ¹³C values with low C₄ input in a normal mixed diet. These results are presented at the site-level to account for local variability in threshold δ¹³C values for C₄ consumption (Fig. 5). The ecozone cluster model was used to create a map of interpolated threshold δ¹³C values for C₄ consumption and hence suggest environmentally adjusted threshold δ¹³C values for C₄ consumption (Fig. 6C). The European background maps used to create Figs. 5 and 6 are vector map data from https://www.naturalearthdata.com/, implemented using the package rnaturalearth in R-Software^132,144.

Materials & correspondence

The corresponding authors are MLCD and MK. The isotopic dataset used in this study is available from the open access repository: Depaermentier, M. L. C. (2025). Isotopic Dataset to: Depaermentier, MLC, Kempf, M, Motuzaitė Matuzevičiūtė, G. “Environmentally adjusted δ¹³C thresholds for accurate detection of C₄ plant consumption in Europe” [Data set]. In Communications Earth & Environment. Zenodo. https://doi.org/10.5281/zenodo.17571650 [ref. ²⁷ in this paper]. The data to reproduce the ecozone clusters is available from this open access repository: Kempf, M. (2025): Related files to: Depaermentier, MLC; Kempf, M; Motuzaitė Matuzevičiūtė, G: Environmentally adjusted δ¹³C thresholds for accurate detection of C₄ plant consumption in Europe (2025) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.15695070 [ref. 147 in this paper]. Climate variables used in this article are freely available from Karger et al. (2017): https://chelsa-climate.org/ (last accessed 19th of June 2025) [ref. ¹⁴³ in this paper]. The Digital Elevation Model (DEM) can be downloaded from the USGS earthexplorer server: https://earthexplorer.usgs.gov/, last accessed 19^th of June 2025.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.