The animals originally selected for CSIA-AA from each site had similar bulk collagen N isotope values close to the mean average for each site (ESM 3; Table S1). Consequently, the results of this study must be interpreted as the likely diet of a particular group of animals from each site, which probably had access to similar plant resources and, in the case of the domestic herbivores, had been managed in similar ways to others in the assemblage.
Exploring the relationships between AA δ15N values in modern reference plants
Nitrogen is incorporated into plants via the GS-GOGAT pathway as Gln and Glu (Miflin and Lea 1977), which was confirmed by 15N-labelling experiments, in which plants show the greatest amount of 15N uptake in these AAs (Yoneyama et al. 2003). Since Glx is central to plant AA metabolism due to providing amino groups to the other AAs (Forde and Lea 2007), it is to be expected that most AA δ15N values would be highly correlated with that of Glx, which is found to be the case for Ala, Asx, Leu, Pro, Val and Phe (Fig. 1). A decreased degree of correlation between an AA and Glx is likely only to result when there are many intermediate steps to biosynthesis or when there are partially self-sustaining recycling mechanisms.
Gly and Ser are involved in one such recycling pathway. The SHMT-SGT pathway uses two molecules of Gly to produce one molecule of Ser. Excess Ser is in turn recycled into molecules of Gly. This pathway produces a deficit of one molecule of Gly if Ser is recycled from excess Gly. The Gly deficit can be balanced by biosynthesis from Glu (Bauwe et al. 2010). This is likely to explain why Gly and Ser can be less correlated with Glx than many of the other AAs (Fig. 1). Previous studies have indicated that increased biosynthesis of lignin via the phenylpropanoid pathway can lead to elevated Phe δ15N values (Kendall et al. 2019). This phenomenon has been suggested to be a distinguishing factor between cereal grains and rachis as well as herbaceous and woody plants, since rachis and woody plants contain more lignin than grains and herbaceous plants (Kendall et al. 2019). In Fig. 1, Phe is less correlated with Glx in herbaceous plants compared with woody plants. Increased demand for Phe in woody plants for lignin production would lead to higher levels of Phe biosynthesis from Glu, thereby increasing correlation with Glx. Interestingly, the consistent correlation of Phe with Glx in plants (Fig. 1) could explain the apparent stability of the Δ15NGlx−Phe value (Steffan et al. 2013) commonly used to quantify trophic level differences (McMahon and McCarthy 2016). The implication is that if the correlation between Phe and Glx weakens, such as in the case with herbaceous plants or cereal grains, the Δ15NGlx−Phe value could become less reliable. Finally, in plants, Thr is biosynthesized from Asp without breaking of a C − N bond. Lack of consistent correlation of Thr with Asx (Fig. 1) could be explained due to conversion of Thr to Gly through threonine aldolase (Joshi et al. 2006), the activity of which may be highly variable between species. Our loading plots of modern reference plant AA δ15N values (Fig. 1) is therefore consistent with known AA metabolic pathways and validates the differences in metabolism between various different plant categories under analysis.
These differences in AA metabolism between the different plant categories means that they are reasonably well-separated in a score plot (Fig. 2a). While a loading plot could not be generated for the legumes, the patterning of their AA δ15N values differs significantly from that of the other plant categories (ESM 1; Table S1, Styring et al. 2014a), allowing them to be separated from the other plants in Fig. 2a.
Assessing the potential and current limitations of using AA δ15N values and FRUITS to reconstruct herbivore diet
Modern reference cattle overlap with herbaceous plants in Fig. 2b as expected in the score plot. However, while the overlap demonstrates a strong relationship between herbaceous plants and consumer diet, FRUITS was utilized to provide a quantitative estimate of percentage contribution from various dietary sources. The known diet of the control cattle from North Wyke Farm Platform was 100% C3 herbaceous plants, which is not what the FRUITS model predicts (Table 1).
Given the evident limitation in applying the FRUITS model, a combined percent contribution of crops (grain, rachis, legume) and non-crops (herbaceous and woody plants) was calculated in Table 1. The score plot (Fig. 2a) shows sufficient separation of crops and non-crops, which justifies this separate categorization. In the case of control cattle, this leads to a predicted 57% contribution from non-crop sources and 43% from crops. A 57% contribution of non-crops to the diet remains significantly lower than the actual contribution of 100%. This error can be attributed to the large standard deviations in the source data (as illustrated by data dispersion in Fig. 2a). Moreover, given the large uncertainties, it is unsurprising that there is insufficient statistical certainty to assign 0% contribution to crops. If some value (such as 10%) must be assigned to a category due to statistical uncertainty, even when the real contribution is 0%, then this error implies overestimation of crop contribution in Table 1, given that there are three categories of crops and only two of non-crops. Since the same source plant values were used when analyzing the unknown data, the predicted diets of control cattle outline the limits to which FRUITS results can be interpreted. Moreover, applying FRUITS to CSIA-AA data necessitates estimating percent contribution of different source AA amino groups to biosynthesized AAs found in bone collagen. Without a greater understanding of AA cycling through the ‘metabolic pool’ (O’Connell 2017), which may vary between species and can be influenced by rumen microbes, it is difficult to achieve more precision. While uncertainty can be reduced with priors in Bayesian models, use of bad priors can have significant deleterious consequences for decreasing model accuracy, while giving the appearance of increasing precision (Cheung and Szpak 2020). Given these limitations and the relative performance of PCA and FRUITS in predicting the diet of control cattle, PCA results should be emphasized over FRUITS.
Using CSIA-AA to reconstruct herbivore diet at three Neolithic archaeological sites
Çatalhöyük
The overlap of domestic and wild herbivores with herbaceous plants in the score plot (Fig. 3a) is consistent with current interpretations of animal management strategies at Çatalhöyük. Current evidence suggests that the Konya plain surrounding Çatalhöyük was a dryland anabranching river system in the Neolithic (Ayala et al. 2017) and strontium isotope analysis of sheep tooth enamel suggests that sheep were likely to have been herded on this plain rather than in the more distant upland region (Bogaard et al. 2014). Aurochs, like their domestic counterpart, are primarily grazers (Gordon and Prins 2008). Therefore, in a region such as the Konya plain, it is reasonable to assume that herbaceous plants would have been the dominant dietary component. Likewise, the C and N isotope values of sheep vary greatly (δ13C range 7.4‰, δ15N range 9.5‰, n = 235), which Pearson et al. (2007, 2021) suggests reflects herding in multiple isotopically distinct environments rather than as a result of feeding on crops used as fodder. Moreover, recovered archaeological goat/sheep dung pellets also reflect high dietary variance with weed species from wetlands (Bolboschoenus glaucus, Carex), common weeds (e.g. Chenopodium album, Polygonum aviculare), and C4 species likely derived from saline marshes (Aeluropus, Crypsis, Sporobolus) (Charles et al. 2014). The results in the score plot (Fig. 3a) corroborate these pre-existing sources of evidence.
Furthermore, the bulk collagen δ13C and δ15N values of ovicaprids (δ13C: −19.4 ± 0.9‰, δ15N: 5.2 ± 0.1‰, n = 5) and aurochs (δ13C: −18.7 ± 0.8‰, δ15N: 9.7 ± 1.4‰, n = 5) overlap with each other, which further supports the interpretation suggested by our modeling that domestic herbivores were herded in a similar environment to their wild counterparts (for a complete list of bulk collagen values see ESM 3; Table S1).
Makriyalos
The domestic and wild herbivores from Makriyalos form two distinct groups in the score plot (Fig. 3b). This might suggest tightly controlled herding of domestic herbivores at Makriyalos, meaning that they had access to plants that were distinct from those available to the red deer. The domestic herbivores overlap mostly with herbaceous plants, while the red deer plot closer to woody plants and rachis. These results support the possibility of herding of domestic herbivores in local saline marshes, as suggested in Vaiglova et al. (2018). The more mixed diets of the two red deer individuals are consistent with known feeding patterns of red deer. Contrary to most deer species, the red deer is not predominantly a browser. Red deer obtain as much as a third of their diet from grasses with the rest made up of concentrate foods which include items such as leaves, shrubs, and seeds (Gebert and Verheyden-Tixier 2001). In regions with sparse woodlands such as the plains surrounding Makriyalos, grasses and herbs can make up a significant portion of diet, so a large herbaceous plant contribution to red deer diet is not unexpected.
Vaihingen an der Enz
There is no distinct separation of domestic and wild herbivores in the score plot for Vaihingen (Fig. 3c), and both overlap with herbaceous plants. The red deer did not overlap with the woody plants as might be expected. Instead, red deer appear to have had a similar diet to the domestic fauna. This interpretation agrees with bulk C and N isotope values. Despite relatively high cereal grain δ15N values (Fraser et al. 2013), cattle δ15N values (6.8 ± 0.2‰, n = 5) are low and not significantly different from red deer (6.2 ± 0.1‰, n = 5), which does not suggest crop foddering. Cattle δ13C values (− 22.7 ± 0.5‰) are also similar to those of red deer (− 23.2 ± 0.3‰) and are noticeably more depleted (average difference approximately 4‰) compared to herbivores from the other sites, which is likely a result of the canopy effect in wooded environments (see ESM 3; Table S1). These results suggest that cattle and red deer fed in similarly wooded environments.
The similarity between wild and domestic herbivores can be explained either by cattle grazing in a similar wooded environment to deer or through some degree of leaf foddering as has been suggested for herbivores in Switzerland (Doppler et al. 2017) and France (Balasse et al. 2012). Deer and cattle investigated in this study have lower δ13C values than sheep (–21.8 ± 1.7‰, n = 4), goats (–21.1 ± 0.3‰, n = 2), and domestic pigs (–21.3 ± 0.3‰, n = 23) from the same site (Fraser et al. 2013). This further supports evidence of some woodland browsing or leaf foddering but from isotopic dietary reconstruction alone, it cannot be determined if lower δ13C values are due to leaf foddering or herding in woodland environments.
Complications associated with C4 plant consumption
While this model can separate out clusters of C3 plants with distinctive metabolisms, a significant limitation remains due to the lack of C4 plants in the source data. When compared with C3 plants, C4 plants have somewhat elevated bulk δ15N values (by ca. 1.5‰) in addition to elevated δ13C values (Hare et al. 1991). Furthermore, Hare et al. (1991) showed that the δ15N values of the AAs Gly and Ser are significantly higher in C4 plants compared to C3 plants (δ15N values were elevated by 3.2‰ and 3.7‰ respectively for Gly and Ser in C4 plants). The reference plants in this study are all C3 and the trophic offset between plant diet and herbivore consumer was determined with C3 plants (Kendall et al. 2017), which given δ15N value differences between C3 and C4 plants, is not necessarily directly applicable to herbivores fed with C4 grasses.
These metabolic differences between C3 and C4 plants are encouraging for future studies; however, at the present, it means that for sites with evidence of C4 plant consumption, the models may not be accurate. At Çatalhöyük, it is possible that some C4 plant contribution could have skewed the FRUITS results, as relatively high δ13C values of both wild cattle and domestic ovicaprid bone collagen indicate a dietary contribution of C4 plants (Pearson et al. 2015) and seeds of C4 grasses have been found in sheep dung pellets (Charles et al. 2014). Likewise, herding in saline marshes is also a possibility for sheep and particularly cattle at Makriyalos. The average bulk δ13C value of − 19.4 ± 0.9‰ (ESM 3; Table S1) for sheep (n = 5) indicates that they could have had a mixed C3/C4 diet while cattle at Makriyalos show elevated δ13C values averaging at − 18.0 ± 2.1‰, with the median at − 16.7‰ for individuals analyzed in the score plot (n = 4) (ESM 3; Table S1). More representative sampling of cattle concludes with similar findings (Vaiglova et al. 2018). Therefore, reference C4 plants and unique saline marsh species should be sampled to further augment the model.