Living off the land: Terrestrial-based diet and dairying in the farming communities of the Neolithic Balkans

The application of biomolecular techniques to archaeological materials from the Balkans is providing valuable new information on the prehistory of the region. This is especially relevant for the study of the neolithisation process in SE Europe, which gradually affected the rest of the continent. Here, to answer questions regarding diet and subsistence practices in early farming societies in the central Balkans, we combine organic residue analyses of archaeological pottery, taxonomic and isotopic study of domestic animal remains and biomolecular analyses of human dental calculus. The results from the analyses of the lipid residues from pottery suggest that milk was processed in ceramic vessels. Dairy products were shown to be part of the subsistence strategies of the earliest Neolithic communities in the region but were of varying importance in different areas of the Balkan. Conversely, milk proteins were not detected within the dental calculus. The molecular and isotopic identification of meat, dairy, plants and beeswax in the pottery lipids also provided insights into the diversity of diet in these early Neolithic communities, mainly based on terrestrial resources. We also present the first compound-specific radiocarbon dates for the region, obtained directly from absorbed organic residues extracted from pottery, identified as dairy lipids.


Introduction
The earliest attempts in the domestication of wild plants such as barley, lentil, einkorn and emmer wheat, and animal species (cattle, sheep and goat), were identified in Southeast a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 direct dating of the occurrence of specific foodstuffs, such as carcass or dairy products [33][34][35]. In addition to the conventional radiocarbon dates obtained from some of the sites (shown S5.1 Fig in S5 File), here, we dated two potsherds containing dairy products, providing, for the first time, a direct date for milk use in the Balkans.

Organic residue analyses of pottery
Lipid analysis and interpretations were performed using established protocols described in detail in earlier publications [36,37] (see also S5 File). Interpretable lipids (in concentrations > 5 μg/g of sherd) were recovered from 55 (26%) of the 213 samples, with the recovery rate from each site being 32%, 21%, 26% and 18% from the sites of Starčevo-Grad, Magareći Mlin, Vrbjanska Čuka and Rutonjina Greda, respectively. This is comparable to recent studies in the region, for example, the overall recovery rate from sites in Neolithic Greece was 23% although rates varied between the Early, Middle and Late Neolithic [20]. Furthermore, the rate was 22% overall at five Early Neolithic sites in the Iron Gates region of the Lower Danube [22] and was 22% overall in 7 sites from the Hungarian plain and the Balkans area [16,20,22,25,38]. The majority of the lipid profiles (n = 47) comprise the free fatty acids, palmitic (C 16 ) and stearic (C 18 ), typical of a degraded animal fat, as the most abundant components (Fig 2a) [23,39]. Other lipid classes were detected, comprising aliphatic lipids including n-alkanes and n-alcohols, which will be discussed further below.
GC-C-IRMS analyses were carried out on these 47 absorbed lipid residues to determine the δ 13 C values of the major fatty acids, C 16:0 and C 18:0 , and ascertain the source of the lipids extracted [24,40]. The δ 13 C values obtained for modern reference animal fats, from animals raised on a pure C 3 diet are grouped within confidence ellipses (±1s), onto which the values from the archaeological pottery have been plotted (Fig 3a-3d).
Ruminant dairy fats are differentiated from ruminant adipose fats when they display Δ 13 C values of less than -3.1 ‰, known as the universal proxy [42,46]. Evidence for the processing of secondary products, such as milk, butter and/or cheese, is found at all sites (total, n = 23), albeit in varying amounts, i.e. at Vrbjanska Čuka (n = 1, 17% of lipid bearing sherds at the site), Starčevo-Grad (n = 12, 48% of lipid bearing sherds at the site), Magareći Mlin (n = 8, 62% of lipid bearing sherds at the site) and Rutonjina Greda (n = 2, 67% of lipid bearing sherds at the site). A further 18 sherds plotted in the ruminant carcass region, Vrbjanska Čuka (n = 3, 50% of lipid bearing sherds at the site), Starčevo-Grad (n = 10, 40% of lipid bearing sherds at the site), Magareći Mlin (n = 4, 31% of lipid bearing sherds at the site) and Rutonjina Greda (n = 1, 33% of lipid bearing sherds at the site). No vessels plot solely within the non-ruminant (pig) region although a number of vessels from three sites, i.e. Vrbjanska Čuka (n = 2, 33% of lipid bearing sherds at the site), Starčevo-Grad (n = 3, 12% of lipid bearing sherds at the site) and Magareći Mlin (n = 1, 8% of lipid bearing sherds at the site), plot between the ruminant and non-ruminant region, suggesting some mixing of these products (Fig 3e-3h).
Aquatic/freshwater resource processing. All fatty acid methyl esters (FAMEs) were analysed by GC-MS in SIM mode to check for the presence of aquatic or freshwater biomarkers, namely ω-(o-alkyl phenyl) alkanoic acids (APAAs) and vicinal dihydroxy acid (DHYAs) which originate from the degradation of poly-and monounsaturated fatty acids found in marine or freshwater fats and oils. These are routinely used to detect marine/freshwater product processing [47][48][49][50][51]. Only one potsherd, Vessel ST35, contained the C 18 and C 20 APAAs, but DHYAs were not identified in any sherds. Hence, although freshwater aquatic products may have been mixed with terrestrial products in this vessel, albeit in low abundances, there is no evidence for sustained processing of freshwater resources at these sites.
Plant processing. Interestingly, a number of lipid profiles from Starčevo-Grad (n = 2, ST73, ST83) and Vrbjanska Čuka (n = 3, VC02, VC09, VC10) differ from those typical of animal fats, comprising sequences of even-numbered long-chain fatty acids (LCFAs), containing C 20 to C 30 carbon atoms, generally dominated by the C 24 (Fig 2b). These LCFAs are strongly indicative either of an origin in leaf or stem epicuticular waxes [52][53][54][55] or, possibly, suberin [56][57][58][59], an aliphatic polyester found in all plants. Although primarily found on the surface of plant leaves, sheaths, stems and fruits, epicuticular waxes are also found associated with other plant organs, i.e. seed oils and coats, flowers, bark and husks [54]. However, these LCFAs are not diagnostic to families of plants and so cannot be used as anything other than a general indicator for plant processing.
Also present in the profiles is a series of odd-over-even long-chain n-alkanes, ranging from C 25 to C 33 , generally dominated by the C 29 n-alkanes (Fig 2b), albeit in low concentrations.

PLOS ONE
Living off the land: terrestrial-based diet and dairying in the farming communities of the Neolithic Balkans Alkanes are also common components of waxes, usually occurring in low concentrations [60], although occasionally they are the dominant lipid, e.g. the leaf wax of Cotyledon orbicularis is almost entirely comprised of alkanes [61]. Long-chain n-alkane distributions occur in the range C 25 to C 35 [62], with an odd-over-even predominance [61]. The dominant chain lengths vary across plant taxonomic groups but the C 27 , C 29 , C 31 and C 33 homologues usually predominate [63]. Significantly, an analysis of leaf wax alkanes extracted from 93 species belonging to five subfamilies, Bambusoideae, Pooideae, Arundinoideae, Chloridoideae and Panicoideae, of the Gramineae (grass family), showed that the C 29 and C 31 n-alkanes dominated [64]. This combination of LCFAs and n-alkanes strongly suggests the processing of plant material, likely leafy plants and/or wild grasses, within these vessels. Interestingly, two sherds from the Vrbjanska Čuka pots (VC09, VC10) comprising plant lipid biomarkers, plot in the ruminant adipose region (Fig 3e, red triangles), suggesting mixing of leafy plants and carcass fats from cattle, sheep or goat, possibly in the form of a stew. In contrast, one vessel from Starčevo-Grad (ST83), which included LCFAs and n-alkanes, plotted within the dairy region (Fig 3f, red  triangle).
Beeswax. A number of lipid profiles (n = 7) comprised series of long-chain even-numbered n-alkanoic acids (C 20 to C 26 ), n-alkanols (C 24 to C 32 ), and n-alkanes (C 25 to C 31 ). These lipid profiles are generally indicative of the presence of beeswax and thus these vessels (four from Starčevo-Grad, ST54, ST55, ST82 and ST90 and three from Magareći Mlin, MM133, MM142 and MM149) were selected for further analysis by solvent extraction [36] to identify higher molecular weight compounds, wax esters, which would confirm the presence of beeswax in the sherds (see S5 File). Of these, the complex mixture of compounds seen across the profiles of ST54, ST55, ST82, ST90, MM142 and MM149 (Fig 4a and 4b) comprises the  The three fields correspond to the P = 0.684 confidence ellipses for animals raised on a strict C 3 diet in Britain [41]. Each data point represents either: blue circle-terrestrial animal product; red triangle-plant/animal product mixture; Plots e, f, g, h, show the Δ 13 C (δ 13 C18:0 -δ 13 C 16:0 ) values from the same potsherds. The ranges shown here represent the mean ± 1 s.d. of the Δ 13 C values for a global database comprising modern reference animal fats from Africa [42], UK (animals raised on a pure C 3 diet) [36] Kazakhstan [43], Switzerland [44] and the Near East [45], published elsewhere.
https://doi.org/10.1371/journal.pone.0237608.g003 following homologous series, C 25 to C 31 carbon number n-alkanes displaying a unimodal distribution possessing a strong odd-over-even predominance and a series of C 24 -C 26 long-chain alcohols in which n-tricontanol (C 30 ) and n-dotricontanol (C 32 ) were the major components. Eluting at longer retention times were a series of C 40 -C 54 carbon number palmitic acid wax esters, confirming the presence of beeswax. A series of hydroxy palmitic acid wax esters, eluting at somewhat longer retention times than the wax esters, in which the C 46 and C 48 homologues were the most abundant components, are also present in samples ST54, ST55, ST82 and ST90. Sample MM133 did not contain any wax esters, simply a series of odd-numbered nalkanes, and therefore cannot be unambiguously confirmed to contain beeswax. In summary, four potsherds were unambiguously confirmed to have contained beeswax, and two others highly probable.
Fresh beeswax comprises a complex mixture of aliphatic compounds consisting of series of homologues differing in chain-length by two methylene groups. Medium-chain nalkanes range from C 23 to C 31 (with C 27 dominating in A. mellifera), and n-alkanoic acids from C 20 to C 36 (usually dominated by lignoceric acid (C 24 )). Monoesters comprise predominantly alkyl palmitates (C 38 to C 52 ), with characteristic hydroxy monoesters comprising long-chain alcohols (C 24 to C 38 ) esterified mainly to hydroxypalmitic acid, ranging between C 40 and C 54 [65]. Although it is relatively resistant to degradation, the chromatographic profile of ancient beeswax often presents significant differences to that of contemporary beeswax [66,67]. For example, the free n-alkanols do not occur in fresh beeswax but are found in aged wax, due to hydrolysis of the wax esters. Furthermore, a preferential loss of lower carbon number n-alkanes may induce a modification of the n-alkane profile through time [68].
Finally, vessel ST54, displaying a typical beeswax profile, also contained high abundances of the major fatty acids, C 16:0 and C 18:0 , typical of animal product processing (Fig 4a). This suggests that beeswax/honey may have been mixed with animal products in this vessel. These were also present in low amounts in vessels ST55 (Fig 4b) and ST82.
Direct dating of dairy residues. Two potsherds from the sites of Starčevo-Grad (ST118) and Vrbjanska Čuka (VC26) containing absorbed lipid residues in high concentrations were radiocarbon dated using a compound-specific approach (CSRA) [33,34]. The method employed is based on the extraction of C 16:0 and C 18:0 fatty acids residues from the clay matrix of potsherds and their isolation from exogenous contamination in the pots. In this case, the direct dating of dairy lipids from the potsherds, identified through their δ 13 C values, allows the date of use to be obtained.

PLOS ONE
Living off the land: terrestrial-based diet and dairying in the farming communities of the Neolithic Balkans

Taxonomic composition of faunal assemblages and cattle mortality profiles
In the three faunal assemblages studied, Magareći Mlin, Starčevo-Grad and Vrbjanska Čuka, the vast majority of remains originated from domestic ruminants-cattle, sheep, and goat (Table 1). Cattle remains were most numerous in the samples from Starčevo-Grad (44.0% of all mammal remains identified to the species/genus level) and Magareći Mlin (60.1% of all mammal remains identified to the species/genus level). This corresponds to the pattern observed on other Early Neolithic sites in the North-Central Balkans and southern fringes of the Pannonian Plain, indicative of a mixed herding strategy which involved ovicaprids and high proportions of cattle [16,27,[70][71][72][73]. Considering the size and body mass of cattle, these animals most likely provided the bulk of protein and were the main suppliers of meat. On the other hand, analysed faunal assemblages from Macedonia [74-76 and references therein] indicate that sheep and goat husbandry was more prevalent in the Southern Balkans. Nevertheless, albeit small, the sample from Vrbjanska Čuka suggests that its inhabitants were herding both cattle (28.7% of all mammal remains identified to the species/genus level) and ovicaprids (50.5% of all mammal remains identified to the species/genus level).
The paucity of the three studied assemblages hinders our understanding of cattle slaughter patterns. At Magareći Mlin and Starčevo, this could be related to the techniques used to recover animal bones (see S5 File). However, in case of Vrbjanska Čuka, both hand-collection and flotation were employed, so the low MNIs could be related to deposition practices, i.e. cleaning of building floors prior to abandonment [77]. Nevertheless, some trends may be discerned. Profiles are constructed for Starčevo-Grad (MNI = 14) and Magareći Mlin (MNI = 15), whereas the minimal number of individuals based on dental elements from Vrbjanska Čuka was insufficient for profile construction (MNI = 4). Nevertheless, given the evidence of milk exploitation and processing at the first two sites by means of lipid residue analysis, it is worth cross-referencing this data with cattle mortality profiles, and, ultimately, completing the assessment using the results obtained by stable isotope analysis.
The cattle from Starčevo-Grad consisted of individuals <6.5 years old (Fig 5a). The number of individuals who had not reached sexual maturity (five infant and six juvenile animals) significantly exceeds the number of adults (three animals). There were five animals aged <1 year, distributed into four different age classes, possibly related to perinatal mortality and both preand post-lactation slaughter. Their presence indicates culling before reaching optimal weight for meat exploitation. The largest number of animals (four), belonging to the age class of 1-2 years, produces a rather significant peak in frequency density and is presumably related to meat exploitation. Among adults, three animals were 2-4 years old, two were 4-6.5 years old, whereas mature and senile animals were completely absent. Given the evidence of ruminant adipose and dairy products in pottery fragments from Starčevo-Grad, the culling profile might be best explained by mixed meat and milk exploitation.
In the cattle assemblage from Magareći Mlin, all age stages except 6.5-9 years are present (Fig 5b). The number of individuals who had not reached sexual maturity (seven infant and juvenile individuals) is close to the number of adults and old adults (eight). There were four animals aged <1 year, each in a different age class, therefore no conclusion could be drawn regarding pre-or post-lactation slaughter. Nevertheless, their presence indicates culling before reaching optimal weight for meat exploitation. Similar to the Starčevo-Grad mortality profile, a significant peak is visible in the age of 1-2 years, presumably related to meat exploitation. The distribution of adults and old adults, probably mainly fertile and lactating females, indicates that stock breeding was not solely oriented towards beef production.

Cattle weaning patterns
Cattle teeth selected for determination of Δ 15 N and Δ 13 C values of dentine collagen originated from individuals aged 4-6 months (LBMM002 and LBSG060) and 12-15 months (LBMM022). The results from the sequential analysis of δ 15 N and δ 13 C values in dentine collagen are shown in S3 File and on Fig 6. Most collagen extracts meet the criterion for good quality preservation (C and N contents, C/N atomic ratio). The sampling performed on the anterior and posterior lobes of each molar shows great consistency in δ 15 N and δ 13 C values (SI 3), suggesting a good preservation of the original stable isotope ratios. The two cattle teeth from Magareći Mlin (LBMM) show decreasing δ 15 N values from the crown apex (earliest formed) towards the root (latest formed). In LBMM022-m2, δ 15 N values around 9 ‰ over the first 25 mm of crown decrease along tooth crown and stabilize around 6.4 ‰ (Fig 6).
The amplitude of the shift in δ 15 N values (3.1 ‰) is similar to the 15 N-enrichment between trophic levels (+3-4 ‰, [78,79] and the +3.2-3.6‰ 15 N-enrichment between cow's milk and diet [80,81]), likely reflecting the weaning process. The δ 15 N values stabilize in the latest 15 mm dentine formed, suggesting that weaning was completed. In LBMM002-m1, δ 15 N values around 8.5 ‰ over the first 10 mm of the molar crown decrease and reach 6.5 ‰ in the dentine latest formed (Fig 6). Although the δ 15 N values have not stabilized and despite lower intra-tooth amplitude of variation (2 ‰), the sharp decrease in δ 15 N values also suggests a well-advanced weaning process. It may be concluded that both these animals were culled after lactation had substantially decreased or ceased. The cattle tooth from Starčevo-Grad (LBSG060-m1) yielded a different pattern with δ 15 N values increasing from 6.4 ‰ in the crown apex to 8.4 ‰ in the latest formed dentine, a value similar to the highest values measured in LBMM002-m1 and LBMM022-m2 (Fig 6). The weaning process is not recorded in this molar, meaning that this calf died or was culled before its mother's lactation ended. In a

PLOS ONE
Living off the land: terrestrial-based diet and dairying in the farming communities of the Neolithic Balkans

PLOS ONE
Living off the land: terrestrial-based diet and dairying in the farming communities of the Neolithic Balkans context where milk exploitation was demonstrated by the presence of milk residues in ceramic potsherds, this result is not in favour of the hypothesis of calves being maintained alive until lactation ends in order to stimulate milk let down [30, 82,83]. However, this isolated find prevents any conclusion on this matter. The pattern of increasing δ 15 N values observed in LBSG060-m1 might reflect the shift from in utero to sucking diet. In modern cattle, the first molar crown is already partly formed at birth [84]. Yet, although a prenatal signal is expected in the m1 crown, it was not detected in previous studies in modern and archaeological cattle aged 12 months or more [28,29]. This would confirm the hypothesis of the earliest stable isotope ratios in dentine collagen being partly to completely overwritten by signals acquired during subsequent tooth development stages [29]. An important consequence would be that the age at death is an additional factor to the initial timing of tooth formation-and potentially the major factor to be consideredwhen interpreting stable isotope ratios in dentine collagen in relation to age.
The pronounced difference in the pattern of variation of δ 15 N values between LBMM002-m1 and LBSG060-m1, assigned to the same age class, is enigmatic. It could refer to a radically different developmental story either in the timing of tooth growth or in suckling behaviour. The limited number of teeth sequentially sampled in the present study, and in all previous studies [28-30] prevents a clear understanding of intra-tooth variations in dentine collagen δ 15 N values.
The dentine collagen δ 13 C values vary between -22 ‰ and -20.7 ‰ (Fig 6). These are related to diet and metabolic processes involved in digestion. In prenatal life, the young calf receives carbon derived from its mother's ruminant digestion. When the calf is born, its own digestive system is similar to a non-ruminant digestive system during the first weeks of life, until the consumption of dry food stimulates the rumen development. These transitions

PLOS ONE
Living off the land: terrestrial-based diet and dairying in the farming communities of the Neolithic Balkans between ruminant/non-ruminant/ruminant digestive systems occur as the first molar develops and are suspected to impact δ 13 C values in enamel bioapatite [85]. However, there is little chance that they would impact collagen δ 13 C values, as a survey of large ruminant versus nonruminant wild herbivores has shown that differences in digestive strategy have no influence on the response of collagen δ 13 C values to diet [86]. Consequently, the variations observed in δ 13 C values along the crown of these cattle m1 and m2 are likely to directly reflect changes in the δ 13 C values of diet. Considering a 5 ‰ 13 C-enrichment between the protein fraction of diet and collagen [87], they refer to plant values around -27 ‰ to -25.7 ‰, typical for C 3 plants in an open environment. Seasonal variations of 1 ‰ to 2 ‰ in the δ 13 C values of C 3 plants have been reported, with the highest values occurring during the summer and the lowest during the winter [88]. Herbivores feeding on these plants inherit these seasonal variations, which are commonly observed in sequential δ 13 C values measured in enamel (see for example [89]). The changes in δ 13 C values in the dentine collagen of the Magareći Mlin and Starčevo-Grad cattle might well reflect seasonal variations. If this assumption is correct, all three cattle recorded a late summer/autumn signal in the last dentine formed, meaning that they died at or shortly after this time of the year.

Dental calculus proteins
Proteins were recovered from all dental calculus samples, although total protein identifications varied substantially across the dataset, ranging from a minimum of 8 to a maximum of 319 (mean = 93.5 median = 78, considering only 'supernatant & pellet' extractions) (see S4.2 Table in S4 and S5 Files). The overall protein recovery in this study was low, when compared to previous dental calculus studies applying a similar protein extraction technique. For example, Hendy et al. [32] identified an average of 400 proteins from the dental calculus of postmedieval skeletons from Britain, while Mays et al. [90] identified an average of 140 proteins from Middle Neolithic and Middle Bronze age skeletons from the site of Stonehenge. The relatively low quantity of recovered proteins is likely not related to the amount of starting material used for analysis, as we detected no correlations between the quantity (mg) of dental calculus analysed and the total number of identified proteins within that sample (Fig 7), consistent with observations by Mackie et al. [91] and Hendy et al. [32].
Microbial proteins accounted for the majority of the identified protein families, ranging from 0-88% of the identified protein families (mean 56%) (Fig 8), consistent with the fact that dental calculus is a microbial biofilm. Mammalian proteins, the majority of which were assigned to the human host, accounted for 5-28% of the identified protein families, similar to the proportion of mammalian proteins observed in Hendy et al. [32].
Any proteins reported as contaminants in previous dental calculus studies [32,92] and/or detected in any extraction or instrument blanks were considered potential laboratory contaminants. Overall, the calculus proteomes in this study demonstrated an unexpectedly high proportion of contaminant proteins, representing 7-83% of the identified protein families (mean 30%). The proportion of putative contaminants was particularly high in samples with low numbers of identified proteins families overall (i.e., samples with fewer than 80-100 protein families) (Fig 9). Proteins derived from human skin (e.g., keratin, collagen, hornerin) were particularly abundant across the dataset, potentially resulting from previous handling of the skeletal remains, and/or a lack of preserved endogenous ancient proteins within the calculus. Although 8 individuals displayed evidence of putative plant and fish dietary proteins (S4.3 Table in S4 File), none of the individuals displayed confident proteomic evidence for milk consumption.

The importance of milk and dairying in the northern Balkans
Within the Balkan Peninsula, it is clear that milk and dairying practices become more commonplace as farming spreads from South to North. On one hand, even though present, dairy fats are less frequent at the Neolithic sites of Greece [20,25], which correlates well with central and eastern Anatolia [25,93]. On the other, pottery containing dairy products are much more common at northern Balkan sites [16,18,38]. Our results fit well in this general pattern of the relative importance of dairying in the region. At the southernmost site of Vrbjanska Čuka, not far from the Macedonian-Greek border, only one of the six pottery vessels that contained animal fats was used for dairying, whereas the remaining five pots were used for animal carcass product processing. In the Pannonian plain, at the Starčevo-Grad site 12 (48%) of the lipidyielding vessels contained dairy lipids, while at Magareći Mlin, 8 (62%) of the lipid-yielding vessels were used for dairy product processing. At Rutonjina Greda, two of the three lipidbearing sherds yielded dairy lipids (Fig 10). The processing of milk clearly became an important activity within the subsistence base of the farming communities of the temperate Neolithic in the Balkans, although possibly varying in intensity between different communities. Dairying has already been suggested by Ethier et al. [16] (see also [17,18,25]) as a strategy of early farmers for exploiting an important resource rich in protein and fats, in a challenging continental environment where the Mediterranean suite of domesticates had to be modified. In addition, some traditional dairying practices significantly decrease the lactose content [94][95][96], and the short-lived milk is converted into a durable and storable product.

PLOS ONE
Living off the land: terrestrial-based diet and dairying in the farming communities of the Neolithic Balkans In a wider geographic perspective however, this pattern should be considered with caution, as there are significant regional exceptions. The most striking is the evidence from NW Anatolia (the Marmara region of Turkey), where the overwhelming majority of the pottery contained dairy fats and there was also evidence of prolonged heat processing of milk [25,97]. Interestingly, these sites also had higher proportions of cattle bones in the faunal assemblage than the surrounding regions. Not as pronounced, but significant percentages of dairy lipids were retrieved from the western Mediterranean as well [38,98]. At the same time, it seems that the first LBK (Linearbandkeramik) farmers of the continental regions of central Europe did not rely significantly on dairying [99,100]. Clearly, there are a number of different factors which contribute to shaping early farming economies. Through the studies so far, a highly varied picture emerges, suggesting a more geographically focused, regional approach is best suited for future research.
Returning to the importance of milk, the metaproteomic analyses of dental calculus samples from the same northern Balkan areas did not provide any evidence of consumption of fresh, unprocessed milk. ß-lactoglobulin, a protein found within the whey fraction of milk, is one of the most common dietary proteins observed within previous analyses of archaeological dental calculus [31,32,101,102]. Because it is exclusively found in milk it is a specific biomarker for this product, and it also appears to be a particularly robust protein and is more resistant to enzymatic degradation and microbial proteolysis than other milk proteins [103].
Neither ß-lactoglobulin nor any other milk proteins could be confidently detected within this dataset. Although skeleton SAJ19 displayed a single peptide matching to bovine ß-lactoglobulin, and skeleton LV17 displayed a single peptide matching to Alpha-S1-casein (S4.3 Table in S4 File), these identifications do not meet the criteria of having two or more distinct,

PLOS ONE
Living off the land: terrestrial-based diet and dairying in the farming communities of the Neolithic Balkans and confidently identified peptides, criteria necessary to reduce the likelihood of misidentifications and false positive results [104][105][106]. Considering that ß-lactoglobulin is the dominant protein in the whey fraction of milk, the lack of ß-lactoglobulin in the dental calculus samples cannot be simply explained with the overall poor preservation of proteins. ß-lactoglobulin appears to be a relatively robust protein and has been identified in other Bronze Age and Neolithic samples, where total protein identifications were similar to those observed in this study. For example, even in relatively poorly preserved early Neolithic dental calculus samples from Britain (with total protein identifications ranging from 15 to 128 proteins), ß-lactoglobulin was detected in six of the 10 tested individuals [107].

PLOS ONE
Living off the land: terrestrial-based diet and dairying in the farming communities of the Neolithic Balkans The lack of direct evidence for milk consumption through our dental calculus protein analyses, however, should not be used in support of a straightforward claim that milk was not used in a raw, liquid form. Further studies are needed on a larger group of samples from the region, as well as an assessment of the degradation processes of the milk proteins in the specific postdepositional and post-excavation conditions. At present, we cannot reject the hypothesis that humans consumed fresh, unprocessed milk (or the whey fraction and related products), based solely on our dental calculus results. Animal milk could have served as a substitute for breastmilk in the diet of the youngest members of the community, as baby weaning food. The emergence of new foodstuffs (i.e. milk and ground cereals) would have provided Neolithic mothers (as well as other members of the community) with novel options for feeding infants and small children prior to introducing solids in their diet. Such baby gruels were possibly served by bone spoons made from cattle metapodial bones, artefacts ubiquitous in the Neolithic Anatolia and the Balkans, as a recent study of bone spoons from the site of Starčevo-Grad has confirmed [108]. These utensils bore numerous traces of use and damage which corresponded to milk teeth marks, likely produced by children during feeding and/or chewing.

Animal husbandry strategies, cattle weaning and possible implications for milk exploitation
Archaeozoological analyses of faunal samples from a number of Early Neolithic sites in the Balkans, including those from Vrbjanska Čuka, Starčevo-Grad and Magareći Mlin discussed in this paper, indicate that animal husbandry (and to a lesser degree hunting, fishing and shellfish collection) represented the main economic activities. Similar to the lipid residue evidence, there seem to have been pronounced regional differences in animal husbandry practices, mainly between the southern parts of the Balkan Peninsula, and its central and northern parts, bordering with the Great Pannonian plain. In the former, herding strategies were mainly

PLOS ONE
Living off the land: terrestrial-based diet and dairying in the farming communities of the Neolithic Balkans oriented towards ovicaprids, whereas in the latter, while ovicaprids remain important, a significantly higher proportion of domestic cattle is documented [16,27]. The greater reliance on cattle in the temperate Balkans Neolithic corresponds to the increased presence of dairy fats in the northern parts of the peninsula, which suggests that dairying most likely involved cow's milk (although milk from other ruminants should not be excluded), which resonates well with the evidence from the Marmara region [25]. Furthermore, cattle mortality profiles from Starčevo-Grad and Magareći Mlin are indicative of culling strategies oriented towards both meat and milk exploitation. Isotopic analyses of cattle dentine collagen indicate that two calf individuals from Magareći Mlin (from age classes 4-6 months and 1-2 years) were kept alive until the end of lactation. If a similar pattern is to be assumed for the herd as a whole, this implies that the herders at Magareći Mlin had to share milk production with the calves. On the other hand, at Starčevo-Grad, a single calf individual (from age class 4-6 months) was not weaned at the time of death. It remains difficult to interpret these discrepancies given the paucity of the sample, but they could be indicative of different culling strategies at the two sites, or perhaps of different treatments of individual animals and different circumstances related to slaughter events.

Subsistence diversity in early farming communities (processing of fish and cereals, and beeswax utilization)
Finding proof for dairying in the Balkans is important because milk is associated with early domestic animals, and therefore with the early farming economy. However, the subsistence base of the Neolithic people was far more diverse. At two of the sites (Vrbjanska Čuka and Starčevo-Grad), in five vessels, we were able to identify compounds indicative of cooking plants. Two of them did not contain traces of animal product lipids, which suggests that they were used exclusively for processing plant-based food.
Beeswax was identified in vessels at sites in the temperate northern zone, near the Pannonian plain (Starčevo-Grad and Magareći Mlin). Although beeswax could only be unambiguously identified in four vessels, it seems that four of the likely six pots displaying beeswax residue were mixed with animal products, likely dairy. The other two contained only beeswax. None of these six vessels contained carcass animal fats. It is possible that beeswax was used as a technical solution in the pottery making process, to reduce vessel permeability, but if that was the case, we would expect to identify beeswax in a much larger portion of the assemblage. It is more likely that we are detecting honey as part of the diet. The size of our sample prevents us from suggesting culinary practices or dietary preferences but combining data from a larger assemblage from a geographically constrained location could be an interesting line of investigation.
In our study, we could not detect a firm evidence for aquatic resources on a molecular level, despite the proximity of three of the four sites to the Danube. Thus, the lipid analyses suggest a diet based mainly on terrestrial resources, which corresponds well with the data gathered for the early Neolithic of the wider region [16,17,22]. So far, the only exception from this pattern in the Balkans are the well-known sites from the Iron Gates (less than 100 km SE from Starčevo-Grad), where aquatic products played dominant role in the human diet. This was confirmed not only by organic residue pottery analyses [22], but also by isotopic analyses on human remains [19,[109][110][111] and archaeozoological studies [112][113][114]. Such a diversity in the human subsistence between the specialized fishing communities in the gorges on one side, and the farming villages in the plains on the other reflects the complex dynamics behind the neolithisation process in this specific area, where the interaction between immigrating farmers and local hunter-gatherer-fishers, and the availability of arable land played crucial roles [19,111,115]. The presence of hunter-gatherer-fishers outside of the gorges is still under investigation.

Conclusions
These combined results from bioarchaeological investigations of pottery and human and animal remains provided new information about the early Neolithic diet in the Balkans. Mortality profiles of domestic animals, established through taxonomic assessment of animal bones, as well as isotopic analyses of cattle teeth targeting the detection of weaning patterns, have suggested a mixed stock-herding strategy, i.e. exploitation of both meat and milk. The exploitation of milk was evidenced directly with the identification of dairy fat molecules preserved in the Neolithic pottery. The authenticity of the lipid residues was further confirmed by direct compound-specific radiocarbon dating of the dairy lipids, confirming their Neolithic origin. The presence of dairy residues in the pottery vessels suggests that dairying was a common practice at the onset of farming in the Balkans. It becomes obvious however, that dairying did not have the same importance for the people living in different climatic and environmental conditions. Even though practiced in the Sub-Mediterranean regions, dairying gained more importance in the colder Continental North. The heat treatment (in pottery vessels) during dairying facilitates the separation of the whey fraction (with the lactose and the soluble proteins) from the fatty fraction and insoluble proteins, such as caseins. The whey was probably not used in adult diet, while the fatty curd was transformed into durable dairy products and consumed. The proteomic analyses of dental calculus from humans could not confirm that raw milk (i.e. fresh, unprocessed, whole milk) was also part of the diet. Due to uncertainties surrounding the taphonomy of dental calculus proteins however, the absence of milk proteins in our samples is insufficient to reject the idea of milk consumption. Our residue analyses strongly suggest milk was part of human subsistence and could quite conceivably have been used during the weaning period of babies and young children.
The exploitation and inclusion of different kinds of domesticates in the diet, namely plants, was also confirmed by molecular traces in the pottery residue. This is not surprising, since other evidence, such as carbonized grains, wheat ear imprints, or sickle blades are common finds at Neolithic sites in the region. Nevertheless, this is further direct, unambiguous evidence that plants, including cereals, were cooked as food. Beeswax (honey) traces were also found preserved within the pottery matrix. Finding virtually no trace of aquatic resources in our residue analyses, we conclude that at least as far as cooking in pottery vessels are concerned, the subsistence at these four sites was firmly based on terrestrial resources.