Blind dates: Exploring uncertainty in the radiocarbon evidence on the emergence of animal husbandry in the Dutch wetlands

.


Introduction
The transition from hunting and gathering to arable farming and animal husbandry has captured immense interest since Gordon Childe popularised the "Neolithic transition" (later termed "revolution") as a major watershed in human history and, hence, a key foundation of human society (Childe, 1925(Childe, , 1952)).Today, the emergence and importance of animal husbandry and crop cultivation are still muchdebated topics, notably with respect to the timing and duration of the transition.Acknowledging that this occurred within different windows of time in various geographical regions, research on the transition in the Netherlands has particularly focussed on the introduction and adoption of domesticated animals in the wetlands by hunter-gatherers.In reference to the influential Neolithization model presented by Zvelebil & Rowley Conwy (Zvelebil and Rowley-Conwy, 1984), the Availability model, several scenarios have been presented with different perspectives on the timing and duration of the process (Louwe Kooijmans et al., 2011;Raemaekers, 2003).
However, these scenarios rely on a handful of archaeological sites and radiocarbon dates, which is particularly troublesome when it comes to the exploration of the initial phases of animal husbandry.While direct dates on bones of sheep/goatspecies not indigenous to north-western Europeare available (Çakırlar et al., 2020), these provide no information on the domestication of aurochs/cattle (Bos) and wild boar/pig (Sus), species which are indigenous and well represented within various bone assemblages and are more likely to have been subject to local domestication processes (Çakırlar et al., 2020;Clason, 1977;Ottoni et al., 2013;Park et al., 2015;Peeters et al., 2015;Prummel et al., 2002;Prummel & Niekus, 2011;Scheu et al., 2008).Next to problems of species identificationi.e.distinguishing wild from domesticated animalsthe lack of reliable site chronologies is problematic.Generally, the chronological framing of stratigraphical sequences is based on small numbers of radiocarbon dates, most often on peat, charcoal and archaeological finds, in addition to sedimentary data and layer succession.Information on sampled materials varies, as does the reliability of the respective contextual information.
Taken together, these limitations make it difficult to evaluate the validity of the various scenarios.In this paper, we lay the foundations for future research by reconsidering the available radiocarbon evidence and the potential of legacy data, in particular.First, we identify and briefly discuss some key sources of uncertainty with respect to the reliability and precision of legacy radiocarbon dates in the studied datasets.Next, we assess the reliability of available radiocarbon dates as a basis for the application of Bayesian statistical techniques, in order to improve chronological models at a site level.This is followed by a case study, presenting three Bayesian models for the site of Hardinxveld-Giessendam De Bruin, a keystone site in the Netherlands in the context of the early Neolithisation debate (Fig. 1).Finally, we evaluate the implications of the modelling output for the originally proposed phasing of this site, and make some suggestions for methodological improvement in the chronological modelling of site stratigraphy.Through this approach, we demonstrate the importance of a critical evaluation of legacy data and the urgent need for further chronological studies of these sites.

Sources of uncertainty
The Swifterbant culture, represented by sites across the Netherlands, the Scheldt valley in Belgium, and north-western Germany, has been in the centre of the debate on Neolithization in northwest Europe, since it features the introduction of pottery, domesticated animals and cereal cultivation in the 5 th millennium BCE (Louwe Kooijmans, 1993;Raemaekers, 2003Raemaekers, , 2014;;Raemaekers et al., 2021).During this period, these subsistence changes were so small scale in a largely hunter gatherers' lifestyle, that a new term was coined -the 'extended broad spectrum economy' (Louwe Kooijmans, 1993).In such find assemblages, it remains challenging to get a clear insight into the timing and scale of food production, as well as its significance in a broader sociocultural context.Unsurprisingly, the debate is still ongoing (Raemaekers et al., 2021), with the current dominance of a Long Neolithization model (Zvelebil, 2001) opposed by Short Early (Raemaekers, 2003) and Short Late models (Raemaekers, 2014).
Notwithstanding a range of issues surrounding the evidence on the introduction of cereal cultivation (Cappers & Raemaekers, 2008;Raemaekers et al., 2021) the aspect of animal husbandry has its own specific set of challenges.Currently, the timing of the earliest domesticated animals in the Dutch wetlands relies on four direct dates on sheep/goat specimens from Hardinxveld-Giessendam De Bruin and Brandwijk-Kerkhof in the Rhine-Meuse delta (Çakırlar et al., 2020) (Fig. 2).Alongside these remains, several specimens of cattle (Bos) and pig (Sus) were identified as possible domesticated animals based on bone size (Oversteegen et al., 2001;Robeerts, 1995).However, a recent Fig. 1.Sites mentioned in the text: 1. Hardinxveld Giessendam de Bruin, 2. Brandwijk-Kerkhof, 3. Hazendonk 1-3, 4. Swifterbant S2, 5. Swifterbant S3, 6. Swifterbant S4. 7. Hoge Vaart-A27.Both paleo-geographical maps show significant changes in the geological environment between 5500 and 3850 BCE.Vos et al., 2020, adjusted.zooarchaeological and stable isotope study showed that there is no sign of animal management from the perspective of paleodiet and mortality profiles at Hardinxveld-Giessendam sites (Brusgaard et al., 2022b).Awaiting the results of an aDNA analysis, their domesticated status cannot, as yet, be confirmed.The remains of exogenous sheep/goat are much more likely to provide evidence of gift exchange with neighbouring farmers since these are relatively rare and occur in a context with imported pottery and tools such as flint and axes (Amkreutz, 2013;Clason, 1977;Louwe Kooijmans et al., 2001;van Gijn et al., 2001).A similar problem is apparent in the Scheldt valley where, besides sporadic presence of sheep/goat bones in the Swifterbant assemblage (Crombé et al., 2020) there is no conclusive evidence of domesticated local taxa or animal husbandry before the last quarter of the 5 th millennium BCE (Brusgaard et al., 2022a).
A further challenge, and central to this paper, is the lack of precision of radiocarbon dates falling on the calibration plateau in the last quarter of the 5 th millennium.All dates from that time window return broad probability distributions (4300-4000 BCE), hampering the solid chronological assignment of samples (Fig. 2).Some sitesnotably the Swifterbant type sites in the province of Flevolandare particularly affected, leaving no clue about the actual duration of human occupation (Schepers & Woltinge, 2020;van der Waals, 1977).
Even without the calibration plateau, radiocarbon dates contain an inherent coarseness when approached individually.Once calibrated, their probability ranges are often too broad to make them relevant to the questions at hand.To overcome this issue, it is necessary to contextualise the radiocarbon results within a secure site chronostratigraphy, allowing for a significant narrowing down of the date ranges when statistical methods, such as Bayesian modelling, are employed.This requires high quality legacy data from sites with well-preserved stratigraphy and targeted dating across the anthropogenic context (Bayliss, 2009;Bronk Ramsey, 2008).
Legacy dates from the study area concern sample materials which are variably prone to erroneous age estimates caused by reservoir effects (Boudin et al., 2010;Cook et al., 2001;Philippsen, 2013); contamination due to preservatives, e.g.'archaeoderm' (Slupik, 2000); inbuilt age (Schiffer, 1986); post-depositional disturbance; as well as inadequate statistical analysis of radiocarbon measurements.Despite everincreasing awareness of these issues, current narratives still rely on chronologies based on such data.Moreover, numbers of radiocarbon dates per site and intra-site coverage vary greatly (see App. B).This diversity in the size of datasets leads to an imbalance in the representation of site chronostratigraphies, which has vast consequences on the debate.
This situation, combined with the lack of firm evidence and results from multi-proxy studies leaves the Netherlands as a conspicuously blank area on the Neolithisation map of Europe and thus evokes an urgent need for further research.To be more precise, investigating animal husbandry requires multiple lines of evidence among which understanding chronology is key.Firstly, to understand the start of this process, precise dates on preferably bones of securely identified earliest domesticated animals are needed to determine the 'first' appearance in the region.Secondly, to reconstruct the timeline of the full integration of animal husbandry into subsistence strategies (e.g. from the earliest appearance to the significant reliance on domesticated animals), direct dates on securely identified domesticated animals have to be combined with zooarchaeological data and multi-proxy research, namely paleogenomics and stable isotope analysis.And lastly, the gathered data should be interpreted by using appropriate statistical methods to maximise the contribution of acquired results, e.g.Bayesian modelling of radiocarbon dates which have to be contextualised within reliable site stratigraphies.

Bayesian modelling and data assessment
The preceding overview demonstrates the limitations of using (legacy) radiocarbon dates to gain insights into the timing of early animal husbandry in the Netherlands.However, Bayesian statistical modelling, in combination with 'chronometric hygiene' (Pettitt et al., 2003;Spriggs, 1989;Waterbolk, 1971) has the potential to contribute significantly to this challenge, bridging the gap between legacy data and future chronological analyses.In this section, we briefly outline the basic principles of Bayesian statistics in chronological modelling, as well as the assessment of data reliability.

Bayesian modelling
Bayesian modelling has become popular in archaeological research (Bayliss, 2009;Bronk Ramsey, 2008;Buck et al., 1992), and is widely used via online radiocarbon calibration platforms such as OxCal Bronk  (Bronk Ramsey, 2009a;Reimer et al., 2020).Ramsey, 2009b) and BCal (Buck et al., 1999).The upsurge in the application of Bayesian statistics is largely the result of the fact radiocarbon dates are not Normally distributed (Gaussian), but represent irregular probability densities that are usually expressed as date ranges to 68 or 95 % probability.Classical statistical approaches are not well suited to such non-parametric data (Bronk Ramsey, 1994, 1997;Buck et al., 1992).Another strength of Bayesian analysis is its ability to incorporate contextual information into the calculations, permitting the development of coherent models that can lead to a narrowing down of date ranges.Bayesian analysis can also help to estimate dates for events that were not even sampled, such as transitions, and the beginnings and ends of cultural phases (Bronk Ramsey, 2009b).Finally, it allows dates to be simulated in order to predict the required number of dates to achieve a desired precision.When properly applied, this latter function leads to the most cost-effective sampling strategy (Bayliss, 2009;Bayliss et al., 2007).
In principle, the contextual information encompasses a 'prior' probability distribution.Prior probabilities convey patterns like the super-positioning of layers, or the intersections of features (e.g.pits) -in other words, stratigraphical phenomena which result from processes/ events at different times (the one must be older/younger than the other).Defining prior information is of crucial importance as it directly influences the model outputs.Hence, a stringent sample selection from available datasets (chronometric hygiene) and a sound understanding of the stratigraphical context is required (Bronk Ramsey, 2009b;Pettitt et al., 2003).

Sample reliability (chronometric hygiene)
Principles of chronometric hygiene have been widely applied on various legacy datasets (Napolitano et al., 2019;Pettitt et al., 2003;Taché & Hart, 2013, etc.).Building on the seminal work of Waterbolk (1971), many grading systems have been developed over the years for the assignment of reliability grades to radiocarbon samples, allowing for a systemic re-evaluation of available radiocarbon datasets.The basic goal is to distinguish reliable and contextually sound dates on shortlived terrestrial material from the remainder, which may have incorporated an offset to their contextual age.
Here, we apply the tenets of chronometric hygiene to six archaeological sites of the Dutch wetlands (App.A and B), which were selected on the basis of two criteria: (1) earliest known presence of (possibly) domesticated species, and (2) well-established archaeological stratigraphy and context of available radiocarbon dates.These sites hold potential for future studies in high-precision chronology.Radiocarbon dates pertaining to the selected sites (App.A and B) were classified into four groups according to their level of reliability of the sample material (Table 1).The taphonomy and the context of these dates, which are of vast importance, are discussed further below.
Table 2 shows the grade distribution for the six datasets.First-grade samples are only present at three sites, while the remaining three solely have second, third and fourth grade samples, regardless of the availability of the short-lived material.Charred seeds are present across the stratigraphic profiles at the Swifterbant S2-S4 sites but these materials were dated only at S4 (Cappers & Raemaekers, 2008;Schepers & Woltinge, 2020;van der Waals, 1977) and S3 (Raemaekers, 2015).Also, one should bear in mind that various types of materials (e.g.cereal grains) might only correspond to specific parts of a site's history, implying that a strict focus on sampling such materials can lead to biased chronologies.The species and the taphonomic context are of paramount importance to the reliability of radiocarbon measurements.Such information would entail, for example, whether the sample came from a bulk (e.g.peat sample) or a single charred grain.However, legacy data are often not as detailed which limits its contribution to the chronological analyses.Moreover, the number of radiocarbon dates per site varies greatly, regardless of the availability of short-lived material (e.g.Table 2).For example, the number of dates from S2 is insufficient for any conclusion regarding the time span of human activity.While containing more radiocarbon dates, the datasets of S3 and S4 remain greatly affected by the calibration plateau.This means that these sites require, besides more reliable dates, a precise stratigraphical context, ideally from vertical sampling across the occupational layers, combined with Bayesian modelling.Similar to S2, both Brandwijk and Hazendonk have no dates obtained on first-grade samples to help anchor the chronologies and understand the potential offset of second and fourth grade samples.Hardinxveld-Giessendam de Bruin, while having the highest number of dates, still has over half of these based on lower-grade samples, notwithstanding the taphonomical and contextual challenges discussed below.Hence, a careful reassessment of the currently established (inter-) site chronologies is needed for all sites.

Case study of Hardinxveld-Giessendam de Bruin
Having the highest number of radiocarbon dates available, almost half of which based on first-grade samples, the site of Hardinxveld-Giessendam De Bruin (hereafter 'De Bruin') also entails the most comprehensive chronostratigraphy among all the sites addressed in this paper, resolved through a combination of spatial analysis, radiocarbon dating and inundation history (Mol & Kooijmans, 2001).Moreover, this site is key to addressing questions related to the initial Neolithisation process, especially in terms of animal husbandry within the Netherlands.Thus, it is the most suitable site to explore the possibilities of applying a Bayesian approach to improve the temporal resolution of key datasets.

Prior information
Within 6 lithostratigraphical units, 28 archaeological layers were identified amid a succession of peat, sand and clay (Fig. 3).Most layers were allocated to a single occupation phase, others were assigned to two phases due to unclarity in the lateral variability in sediment facies.This is the case with colluvium layers 40 and 50 which are intercalating into the first and second phase.Some layers expand into several phases, such as layer 32, placed in a hiatus between and partly extending to phase 2 and 3.According to a revised chronology, this hiatus lasts for at least 400 years (Mol & Zijverden, 2007).This significant time gap is followed by a final occupational phase which indicates certain changes in subsistence strategies, prompted by the introduction of the first domesticates and, perhaps, dairy production (Çakırlar et al., 2020;Demirci et al., 2021).

Specifications of models
Following the associated archaeological interpretation, (Fig. 3) we constructed Bayesian chronological models using the program OxCal and the IntCal20 calibration curve (Bronk Ramsey, 2009b; Reimer et al., Material with a (freshwater) reservoir effect Bones of omnivores, aquatic animals and species feeding on aquatic sources, charred food crust on pottery Fourth grade Material with protective coating ('archeoderm') Bones with protective coating which cannot be fully removed 2020) in which we allocated the published radiocarbon dates to their respective occupational phases.These phases are clearly demarcated by the archaeological hiatuses between them.Each Phase2 in OxCal corresponds to an occupational phase (henceforth De Bruin 1, 2 etc.) at De Bruin.These were modelled as a series of sequential Phases enclosed by the Sequence command.Because of relative uncertainty in the succession of layers within the occupation phases at the site, we made the decision to avoid enforcing an order between radiocarbon dates within Phases.This allows for the MCMC process to estimate their most probable age, as well as the dates for the transitions (OxCal Boundaries) between them.Where the context of a radiocarbon sample was ambiguous, it was allocated to a single Phase within a separate Sequence whose Boundaries were cross-referenced to the start and end Boundaries of the main Sequence (e.g.model 3A, Fig. 6).An Outlier probability of 5 % was assigned to the first-grade samples (Bronk Ramsey, 2009b).This feature allows the program iteratively to down weigh results which do not agree with the model as a whole.Inbuilt age (IA) samples, such as wood and charcoal, were treated as likely being older than their archaeological context.Here, successful selections can almost only be made towards younger ages, relative to the unmodelled date.To achieve this, the Charcoal Plus model was implemented.For samples likely to be affected by a reservoir effect, we used the After command which treats them as termini post quem and also enforces a shift to younger ages (Dee & Ramsey, 2014) (Table 3).Finally, the duration of Phases were presented with the Date function, conveniently estimating ranges of Phases.Sum and Interval produced comparable results (see Bronk Ramsey, 2017;Loftus et al., 2016); for Interval and Sum, use the code in Supplementary, App.C).
Models in which Outlier Analysis was not employed were also run as comparisons and sensitivity tests (See App.C and D).Here, we observed the Agreement of the model as an indication of the integrity of the chronostratigraphy.In each case, the result with the lowest individual Agreement was manually removed until an Overall Agreement of 60 % was achieved.When OxCal's Outlier Analysis is employed, however, the Overall Agreement of the model becomes a less reliable measure of the suitability of the overall configuration.In general, the two approaches produced comparable results (see specifications in App.C and D).

Table 2
Number of radiocarbon dates (from archaeological layers and samples directly related to anthropogenic events) per site according to the grade of reliability of the sample material.A detailed list of radiocarbon dates can be found in App.B.

Case study model results
The following models of De Bruin present a time range of over 1000 years which includes periods of human occupation, of low intensity of use, and of hiatuses.Most of the radiocarbon dates are directly associated with archaeological layers while four are from soil monolith tins collected from profile sections, marked as DB in the sample name (see App. B).Here, we present the results of the three main models (1A, 2A, 3A) and refer to the sensitivity test models which are listed in the supplementary information (App.C and D).

Model 1A
Only 11 samples qualified for Model 1A, which constitutes a small sample size considering the timespan of the occupation of the site.While most dates agree with their contextual attribution, we see potential uncertainty in the chronostratigraphy.The plot of modelled probability densities shows an unusually elongated De Bruin 2, comprising two clusters of dates and a possible hiatus in between (Fig. 4).Since this hiatus is not attested archaeologically, we suspect the possible intrusion of younger dates into the end of De Bruin 2. This younger cluster of three dates (GrA-14864, GrA-13278, GrA-62951) is suspicious for several reasons.Sample GrA-62951 was labelled as intrusive in the archaeological report.It was measured on a bone of goat/sheep, and only two bones from such species were found in De Bruin 2, while all the rest were found in De Bruin 3. Taking this evidence in combination with the lack of clear distinction between the layers in this part of the stratigraphy, this attribution to De Bruin 2 was regarded as unsafe at the time of excavation (Oversteegen et al., 2001).Nonetheless, both botanical macroremains samples GrA-13278 and GrA-14864 are allocated to the end of De Bruin 2, layers 32 and 20, respectively.According to the outcome of Model 1A, there was no hiatus between de Bruin 2 and 3, yet there might have been one within the second phase.However, as the latter was not attested archaeologically and yet a longer hiatus between De Bruin 2 and 3 was, the two macroremains dates from the end of De Bruin 2 might also be intrusive material from upper layers.
Our suspicion over the context attribution of these samples derives from the following.Sample GrA-14864 is a peat bulk sample from a monolith tin, containing uncharred Corylus shell fragments.Thus, there is a risk that these might contain a mixture of freshly deposited and intrusive material (e.g.due to rooting), alongside the uncertainty in their anthropogenic status.Furthermore, sample GrA-13278, charred Alnus seeds, was collected from the 'top of phase 2 ′ , layer 32 as a terminus ante quem (Bakels et al., 2001;Mol and Zijverden, 2007).This layer as mentioned above, expands over a hiatus and both De Bruin 2 and 3. Thus, the possibility that botanical sample GrA-13278 corresponds to a hiatus or the beginning of De Bruin 3 is realistic.
Further evidence which suggests context misassociation is provided by botanical sample GrA-13272 (beginning of De Bruin 3), which was deemed completely incompatible with the other dates and the stratigraphic information in both the Outlier-and Agreement-based models.(see Fig. 4) This is because the radiocarbon measurements on this sample suggest a much older date, overlapping with De Bruin 2 and preceding the two botanical remains from lower layers.
In sum, we suggest that there is an issue with the contextual attributions between the end of De Bruin 2, a hiatus, and the beginning of De Bruin 3. Despite this, OxCal could not eliminate the potentially intrusive dates in De Bruin 2. What resulted was an elongation of De Bruin 2, partly because the model was not solidly constrained at the end, since it is followed by a hiatus.Contextually, this lack of clarity between the last two phases at De Bruin is to be expected because the complexity of the stratigraphy intensifies towards the top of the dune, where remains of the last occupational phases are located, resulting in unclear boundaries between archaeological layers.This is a consequence of erosion and multiple reuse of surfaces which occurred at the top of the dune, disturbing the primary archaeological deposits and contextual association, and affecting the end of De Bruin 2 and entire De Bruin 3 (Mol & Kooijmans, 2001).In contrast, the probability ranges for De Bruin 1 seem unchanged compared to the unmodelled dates.However, this resemblance is unsurprising because the sample size here is too small to present a convincing chronology.Therefore, more dates are needed to balance the number of occupational phases with respect to radiocarbon dates.

Model 2A
In Model 2A, we introduced the second grade samples which Fig. 6.Model 3A.Plot with probability density distributions of modelled dates; first, second and third grade samples; Hardinxveld-Giessendam de Bruin.Blue probability densities are estimated time ranges of phases, under Date function (for Sum and Interval, refer to the code in Supl).Light grey present unmodelled while dark grey show modelled probability densities.Samples with reservoir effect are labelled under "Human burial" and "Foodcrust".In the separate sequence, labelled ''Unknown context'', we placed all samples with ambiguous context attribution.
increased the sample size by four new radiocarbon dates, contributing significantly to De Bruin 1 (Fig. 5).The second grade samples, mainly wood and charcoal, seem to fit well with the first grade samples which may suggest that the old wood effect is not so distinguishable in this dataset.Depending on the source of the charred wood (e.g.trunk as opposed to branches and twigs), this offset may, in some cases, be negligible.Some studies suggest that temperate regions have a lower risk of offsets caused by the old wood effect compared to the desert and arid environments, due to the faster decay of organic material which may limit the (re)use of old wood (Kim et al., 2019;Mcfadgen, 1982).Similar indications are evident from dendrochronological studies at the early 5 th millennium site of Hoge Vaart.These studies show a significant decrease in the lifespan of trees between 4800 and 4535 BCE due to the increasingly wet conditions (Peeters, 2007;Vos, 2015).The expansion of alder trees in the area of Hardinxveld sites reflects similar environmental conditions (Bakels et al., 2001).Thus, the offset of charcoal and wood samples in the studied assemblage may have been limited to a few decades at most.However, this is uncertain because of the limited information in legacy datasets, the unknown amount of inbuilt age, as well as its inherent complexity.Hence, to explore this, we need more short-lived dates in all phases, to anchor the material with inbuilt age.
With more radiocarbon dates introduced, the modelled probability densities reached a noticeably tighter date ranges compared to Model 1A (compare de Bruin 1 in Figs. 4 and 5).However, the issue of the two separate clusters of dates in De Bruin 2 was not resolved since no extra dates were introduced within this phase.

Model 3A
In Model 3A we introduced the third grade samples.It is apparent that some of these samples scatter towards ages older than their allocated context, which likely indicates the presence of a freshwater reservoir effect (Fig. 6).This is most pronounced with the two samples of 'foodcrust' (GrA-13317, GrA-13320) in De Bruin 3 which both overlap with De Bruin 2. In the sensitivity test for model 3A (Models 3B, 3C, 3E, 3F, see App.E and F) these samples were labelled as outliers, together with the GrA-13272, suggesting the presence of an age offset.This result is in excellent agreement with recent work on the chemical composition of these same food crusts (Demirci et al., 2021).In this study, lipids extracted from the food remains were found to be associable with freshwater fish, and so it is unsurprising that these dates are likely to be affected by a freshwater reservoir effect.This offset is unpredictable and non-quantifiable, especially in the context of hunter gatherer fisher Mesolithic communities within this area (Boudin et al., 2010).Therefore, we argue that by introducing third grade samples, we raise the risk of introducing erroneous dates, leading to an erroneous model.

Phasing reconsidered
The results of the modelling exercise in this paper suggest a significant chronological offset in archaeological phases compared to previous studies (Fig. 7).In the original publication (Mol & Kooijmans, 2001, pp. 70-71), the duration of phases was inferred from the 95 % (2σ) probability ranges of unmodelled dates, despite the significant overlap between dates from different phases (Fig. 7).Moreover, while some radiocarbon dates were deemed unreliable due to the reservoir effect (e. g. the human burials), samples of charred food crust were treated as reliable.This shifted the related phases in that study to older (likely erroneous) ages.On the contrary, short-lived samples of plant macroremains, which returned dates outside of their expected time ranges, were excluded, such as GrA-13272 (13249 in Fig. 7) and GrA-10950 (13254 in Fig. 7).This led De Bruin 3 to appear several centuries older than in a later interpretation (Table 4, compare the interpretation from 2001 and 2007).Hence, once analysed, the dates on goat/sheep in De Bruin 3 turned up much younger than the time range expected (Çakırlar et al., 2020)  Zijverden , (2007).This date (GrA-64342) coincides well with sample GrA-10950 (13254), which was previously excluded from the time range of phase 3 (Mol & Kooijmans, 2001) (Fig. 7).Short-lived samples from De Bruin 2 which seemed to overlap with De Bruin 3 were also excluded (GrA-14864 and GrA-13278 in Fig. 7).
A similar approach was employed in a chronological model published in 2007 (Mol & Zijverden, 2007, p. 97, Fig. 6); this paper Fig. 8), where human burials were excluded, but charred food remains, likely to contain residues derived from fish and hence susceptible to reservoir effects, were retained as reliable samples.Moreover, radiocarbon dates sampled from the beginning and end of occupational phases were placed outside of the Boundaries of these Phases in OxCal (Fig. 8), where they ought to have been placed inside.While the reason for this might be that these dates were samples at the beginning and end of occupation layers at De Bruin, such allocations of archaeological samples are erroneous in Oxcal and, therefore, resulted in an artificial shortening of the respective Phases.
Once we had accurately translated the archaeological stratigraphy into an OxCal model and included all the dates which have no risk of a reservoir effect (e.g.model 2A), differences in the outputs compared with the model from 2007 were immediately apparent (Mol & Zijverden, 2007, p. 97, Fig. 6) (This paper, Table 4).Moreover, while the absolute chronology of this site was considered to be resolved in some former publications, our re-evaluation indicates that further work is needed, partly due to the improved precision of available facilities (AMS) and partly because the legacy data are limited in its detail.We suggest that future studies should concentrate on the radiocarbon dating of shortlived materials directly associated with secure archaeological contexts.If possible, this dating should be targeted across all occupational phases, with an emphasis towards the of end De Bruin 2 (e.g.Layers 20 and 32) and modelled with solely reliable dates from clear context.
The outcome of the newly proposed phasing of this site affect the perceived timing of the introduction of domesticated animals in the Swifterbant community.So far, 4500 BCE has been identified as the oldest occurrence of sheep/goat species and, possibly, animal husbandry in the Dutch wetlands.Given the stratigraphical ambiguity of some of the dated samples, our findings indicate that the more plausible time period is between 4400 and 4000 BCE, which is significantly younger than the introduction of the sheep/goat in the Swifterbant culture in Scheldt valley (Crombé et al., 2020).However, without Bayesian modelling and strategic dating across stratigraphy, it remains difficult to compare the chronologies of these key sites in the region.

Conclusion
Reconstruction of the timing and duration of past processes, such as the transition to animal husbandry, depends on a firm understanding of site chronologies.Thus, before delving into broader sociocultural questions, it is important to lay a foundation by analysing all available chronological data.In this paper, we presented an overview of the most common issues in the radiocarbon datasets of selected sites related to the early animal husbandry in the Netherlands.We draw attention to the flawed nature of current chronologies, and the underrepresentation of high-quality sample material in the radiocarbon datasets.Issues pertaining the inbuilt age, the freshwater reservoir effect, ambiguous context attribution, and inadequate statistical approaches to radiocarbon measurements impact on the established historical narratives today.Our analysis reveals the urgent need for a revision of existent chronological information.By using Hardinxveld-Giessendam de Bruin as a case study, we have demonstrated the contribution of Bayesian modelling to detecting misassociation and inconsistencies in current chronostratigraphies.Our analysis suggests that the last occupation of this site is around 200 years younger compared to the latest chronostratigraphical interpretation, therefore affecting the debate on the initial stages of Neolithisation of the Dutch wetlands.However, we conclude that this site, being one of the most carefully radiocarbon dated site related to early husbandry in the Netherlands, does not currently meet the required criteria for a reliable site chronology and needs further research.The challenges we have encountered in this case study likely extend to the other sites discussed in this paper, thus endangering the integrity of the current narrative on the chronology of early animal husbandry in the Dutch wetlands.We, therefore, strongly advise reassessment of the legacy datasets as they offer an important foundation for future studies in chronologies.
Funding.This study is part of the project 'The Emergence of Domesticated Animals in the Netherlands', funded by the Dutch

Fig. 2 .
Fig. 2. Probability density distributions of calibrated radiocarbon dates on earliest domesticated species (sheep/goat) in the Dutch wetlands (Çakırlar et al. 2020) plotted against the 2020 calibration curve.Purple distributions on the left fall on a calibration plateau (Bronk Ramsey, 2009a; Reimer et al., 2020).

Fig. 3 .
Fig. 3. Schematic overview of the stratigraphy of De Bruin.The right horizontal pane on the right marks three occupational phases encompassing respective archaeological layers.Figure from Mol and Kooijmans (2001), p.62, fig.3.3.;translated and adjusted for clarity.

Fig. 4 .
Fig. 4. Model 1A.Plot with probability densities of modelled dates; first grade samples; Hardinxveld-Giessendam de Bruin.Blue probability densities are estimated time ranges of phases, under Date function (for Sum and Interval, refer to the code in Supl).Light grey are unmodelled while dark grey show modelled probability densities.

Fig. 5 .
Fig. 5. Model 2A.Plot with probability density distributions of modelled dates; first and second grade samples; Hardinxveld-Giessendam de Bruin.Blue probability densities are estimated time ranges of phases, under Date function (for Sum and Interval, refer to the code in Supl).Light grey are unmodelled while dark grey show modelled probability densities.
Fig. 7. Density distributions of radiocarbon dates (95.4% range) within allocated phases.Black are dates from 2001.Purple represents the new goat/sheep dates (Çakırlar et al., 2020) which align outside of expected time ranges of de Bruin 2 and 3. Notice the significant overlap of dates between phases and hiatuses.Figure from Mol and Kooijmans, 2001, fig.3.5.p.69., translated and adjusted.

Fig. 8 .
Fig. 8. Chronological model of Hardinxveld Giessendam de Bruin according to the 2007 publication, figure from Mol and Zijverden, 2007, Fig. 6, p. 97; redrawn and adjusted for clarity.Probability densities enclosed in red squares are dates placed outside of Phases, contrary to the excavation report where these samples are associated with the anthropogenic occupation (Mol and Kooijmans, 2001).

Table 1
Distinction of radiocarbon dates based on the sample material.

Table 3
Specifications of the 3 main models, with attributed Outlier functions.(for Sensitivity test models, see App.C, Table.C.1.).

Table 4
Comparison of time ranges of phases according to different publications.Time range of phase 2 in brackets are without the three intrusive dates.The last phase (De Bruin 3) is significantly younger, shifting the time range of the introduction of earliest domesticates in the Dutch wetlands.