Early Neolithic population dynamics in the Eastern Balkans and the Great Hungarian Plain

pop-ABSTRACT – In this study, we reconstruct population dynamics in the Early Neolithic of the Eastern Balkans and the Great Hungarian Plain using frequency of radiocarbon dates as a population proxy. The method of summed calibrated radiocarbon probability distributions is applied to a set of dates recently published in Bulgaria and Hungary. The aim is to test the hypothesis of the Neolithic demographic transition (NDT) in these regions and to compare the patterns between these two and neighbouring regions. The results show that episodes of population growth occurred in both regions, which is in partial agreement with the predictions of the NDT theory. Population growth is detected, but it is followed by a bust, rather than stabilisation as predicted for the second phase of the NDT. ra-diokarbonskih (NDT) ki predvideva stabilizacijo.


Introduction
It is commonly accepted that the Neolithic was introduced to Europe from the Near East, and new genetic and bioarchaeological evidence undoubtedly show that these processes included movements of people (Davison et al. 2007;Haak et al. 2010;Brami, Heyd 2011;Fort 2012;Pinhasi et al. 2012;Bori≤, Price 2013;Gurova, Bonsall 2014;Özdogan 2014;Mathieson et al. 2015;Szécsényi-Nagy et al. 2015;Hofmanová et al. 2016).The directions and rates of spread of the Neolithic in different regions of Europe have been important issues in numerous studies, many of which offered possible models of pop-ulation expansion.Different models have been suggested for the process, the most frequently considered being: the wave-of-advance model, leap-frog colonisation, and diffusion of cultural novelties (cultural transmission) (Ammerman, Cavalli-Sforza 1971, 1973;Tringham 2000;Whittle et al. 2002;Bar-Yosef 2004;Pinhasi et al. 2005;Davison et al. 2007;Bocquet-Appel et al. 2009;Bori≤, Price 2013).In recent years, the process of neolithisation has been studied as a more complex combination of demic and cultural diffusion (Wirtz, Lemmen 2003;Fort 2012;2015).
The Neolithic way of life brought changes in subsistence and mobility patterns, and a major shift in population structure and dynamics known as the Neolithic Demographic Transition (NDT).According to Jean-Pierre Bocquet-Appel (2002;2008;2011a;2011b;2013;Bocquet-Appel, Bar-Yosef 2008a), the NDT was a two-stage process -the first stage being characterised by exponential population growth caused by increased fertility, followed by a second stage marked by increased mortality and decelerating growth.It is assumed that the increase in fertility was caused by changes in lifestyle which accompanied the Neolithic -dietary (introduction of new, more nutritious food), and changes in residential mobility (sedentary lifestyle).The increase in mortality that followed, especially among infants, was a result of numerous factors -introduction of new pathogens, lack of drinking water, contamination by feces, reduced breastfeeding and higher workload (Bocquet-Appel 2008.49;2013.2).When the mortality rate equaled the birth rate, population growth stopped.
The pattern of Neolithic population dynamics, based on the frequency of radiocarbon dates from Western, Northern and Central Europe, consists of a population boom that occurred at the beginning of the period, followed by a population bust after a few centuries (Shennan et al. 2013;Timpson et al. 2014).This pattern is in agreement with predictions of the NDT theory regarding its first phase, when population growth should occur.A recent study has shown that a similar pattern is observed in the Early Neolithic of the Central Balkans (Por≠i≤ et al. 2016).
In this paper, we extend the paleodemographic research to the area of Eastern Balkans and the Great Hungarian Plain.The aim is to reconstruct population dynamics and to test the NDT hypothesis on two datasets from Southeast Europe, where such studies have been lacking.
Two routes of expansion of the Neolithic way of life into Europe have been proposed.The continental route led from Thessaly, where Neolithic was introduced around the middle of the 7 th millennium cal BC (Perlès et al. 2013) through the Balkans, and to Central, Western and Eastern Europe (Bocquet-Appel et al. 2009;Brami, Heyd 2011;Özdogan 2014).The maritime (Mediterranean) route led from the coast of the Ionian Sea, along the eastern and western Adriatic coast, and further to the western Mediterranean and Iberia.Recent studies (Wirtz, Lemmen 2003;Bocquet-Appel et al. 2009;Lemmen et al. 2011;Silva, Steele 2014;Weninger et al. 2014;Brami, Zanotti 2015) emphasised the complexity of these processes, and different rates and timings for different European regions.These differences depended on various factors, such as geography, climate and sociocultural trends.
The present state of research suggests that the Neolithic spread to the territory of Eastern Balkans (mostly modern-day Bulgaria), along the valleys of the Vardar, Struma and Marica rivers (Boyadzhiev 2009).After settling in Southwestern Bulgaria and Thrace, Early Neolithic populations spread further north.The Neolithisation of this territory was gradual, in a south to north direction, and absolute dates from sites such as Poljanica-Platoto and D∫uljunica-Sma ˘r-de∏ date the beginning of these processes to around 6200/6100 cal BC (Boyadzhiev 2009.11;Krauß et al. 2014.63, Tab. 1).The Early Neolithic lasted until around 5400/5350 cal BC, with the transition to the Late Neolithic occurring between 5500/5450 and 5400/5350 cal BC, which is thought to have been a gradual and smooth process (Gatsov, Boyadzhiev 2009.26).
The Early Neolithic on the Great Hungarian Plain (modern-day Hungary) is represented by the Star≠evo and Körös cultures.The settling of these early farming populations was concentrated in the valleys of several rivers: Tisza, Körös, Maros and Beretyó.The highest density of sites has been found on the southern part of the Great Hungarian Plain, which is explained as a result of its rich and diverse landscape (Paluch 2012.49).The earliest 14 C dates from the southernmost sites, such as Deszk-1, Pitvaros, Maroslele-Pana and Olajkút, indicate that the beginning of Körös culture can be dated to around 6000/5910 cal BC.A south-to-north expansion is also suggested, with the earliest dates for the northernmost (Upper Tisza) region covering the time span between 5630 and 5470 cal BC (Domboróczki 2010;Domboróczki, Raczky 2010;Oross, Siklósi 2012).

Data and method
In this study, published data from Hungary (Anders, Siklósi 2012) and Bulgaria (Gatsov, Boyadzhiev 2009;Krauß et al. 2014) were used.A total of 179 published radiocarbon dates from 16 Bulgarian sites and 117 dates from 24 sites from the territory of Hungary were analysed (Fig. 1; Appendix 1).Dates with large standard errors (170 radiocarbon years and greater) and dates that were out of the currently accepted chronological range for the Early Neolithic in Southeast Europe ( ~6200-5300 cal BC) were excluded from the analysis.1 1 The population dynamics were reconstructed by applying the summed calibrated radiocarbon probability distributions (SCPD) method.The main assumption of this method is that the quantity of material culture is directly proportional to population size in a certain time interval in a given region (Rick 1987;Shennan et al. 2013;Williams 2012).If the number of radiocarbon dates is large enough, then the frequency of dates from a specific time period will be directly proportional to the quantity of archaeological remains from that period, and hence to the size of the population that produced them.The method that was applied in this study was developed by Stephen Shennan et al. (2013) and Adrian Timpson et al. (2014), and it accounts for biases that can greatly affect the final result: the effects of the calibration curve, research bias and effects of taphonomy.The analysis was performed in R programming language (R Core Team 2014), using the Bchron package for calibrating dates (Parnell 2014), and the INTCAL 13 calibration curve (Reimer et al. 2013).
The research bias is the result of different sampling strategies, depending on particular research questions.In other words, samples are usually not collected randomly, but in order to provide chronological information for specific archaeological contexts.In order to reduce the research bias, a binning procedure was performed at the beginning of the analysis.Radiocarbon dates were binned into site-phases, and sorted in decreasing order within each sitephase.Subdivision into bins within site-phases was performed if the difference between two adjacent dates was greater than 200 radiocarbon years.After the calibration, dates were summed within and between bins, and normalised to produce the final SCPD curve.This procedure controls for research bias due to the different number of dates from different sites and site phases, but it cannot control for the bias resulting from the selection of sites and site phases themselves.The binning procedure performed on the 179 Early Neolithic dates from Bulgaria produced 22 bins, while the same procedure performed on 117 Early Neolithic dates from Hungary produced 26.
Taphonomic bias refers to the loss of archaeological material over time due to various taphonomic factors.In order to address this source of bias, the taphonomic exponential curve equation developed by Todd A. Surovell et al. (2009) was used as a null model.The null model assumes that the population was stationary and that, apart from the shape of the calibration curve, the taphonomy is the only factor which affects the shape of the empirical SCPD curve.According to the probabilities given by the null model, calendar dates from the specified time interval were randomly sampled, which produced a large number of simulated radiocarbon datasets.The number of dates for each simulated dataset is equal to the number of bins in the empirical data set.This procedure was repeated many times; for the Early Neolithic dates from Hungary and Bulgaria, we simulated 10 000 null model SCPDs.Sampled calendar dates were then 'back calibrated' and recalibrated afterwards, and summed to produce the simulated SCPD pattern.Finally, the empirical SCPD curve was compared to the 95% confidence intervals calculated from the simulated SCPD values.When the empirical SCPD is above or below the 95% confidence intervals, there is a statistically significant growth or decline of population relative to the null model.This whole procedure was undertaken in order to assess the statistical significance of the empirical SCPD pattern.

Results
The results of the SCPD method for the territory of Bulgaria are shown in Figure 2.After ~6000 cal BC, the curve begins to increase, reaching a peak around 5700 cal BC, after which it decreases.However, none of the changes in the curve reach the threshold of statistical significance, as the curve is always within the 95% CI limits, meaning that it is consistent with the null model, which assumes uniform population (with effects of taphonomy).
The results of the SCPD method for the territory of Hungary are shown in Figure 3.The empirical curve increases after ~6200 cal BC and goes beyond the upper 95% CI limit between ~5750 and ~5500 cal BC, meaning that in this interval, there was a significant increase relative to the null model, which assumes uniform population (with effects of taphonomy).After this interval, the curve abruptly drops, but stays within the 95% CI limits.Around 5200 cal BC, a minor, but statistically significant drop can be observed.

Discussion
In the Eastern Balkans, the curve started to increase with the beginning of the Neolithic.The peak around 5700 cal BC occurred afterwards, and may be considered as indicative of the NDT.The lack of statistical significance (at the 0.05 level) of the deviation from the null model is most probably due to the low effective sample size (the number of bins is only 22), which implies low statistical power.
The results for the Hungarian Plain dates show a significant peak around ~5750 cal BC, which can be interpreted as population growth at the beginning of the Neolithic, and the signal of the NDT.It is followed by a sharp decrease in the curve at ~5500 cal BC, suggesting a population bust in this period.
The observed pattern -a population increase at the beginning of the Neolithic, followed by a population decrease after about 150 years for Bulgaria, and about 250 years for Hungary -correspond to the boom and bust pattern observed in other regions in Europe and the Balkans (Shennan, Edinborough 2007;Shennan et al. 2013;Timpson et al. 2014;Por≠i≤ et al. 2016;Pilaar Birch, Vander Linden 2017).
When compared to the results obtained for the territory of Serbia (Por≠i≤ et al. 2016), it can be seen that they are quite similar in general, but some regional differences should be further discussed (Fig. 4).In Por≠i≤ et al. (2016) it was shown that the SCPD curve for data from Serbia had two statistically significant peaks ( ~6000 cal BC and ~5650 cal BC) and a major drop between.Two explanations have been proposed.The first perceives these changes as real demographic patterns that reflect major population growth followed by increased mortality or migration, with a rebound occurring after 350 years.The other explanation would be that this result is a consequence of a research bias that led to oversampling the earliest Early Neolithic contexts and that the peak around ~5650 cal BC is most probably the signal of the NDT (Por≠i≤ et al. 2016.6-7).In a recent study by Suzanne Pilaar Birch and Marc Vander Linden (2017), which primarily deals with the correlation between environmental changes and population dynamics dur-

Fig. 3. Results of the SCPD analysis based on the Early Neolithic radiocarbon dates from Hungary. SCPD empirical curve (black line) for Early Neolithic dates, with 95% confidence intervals (shaded) based on 10 000 simulations from the null model (grey dashed line) and 200 year rolling mean (red line); number of dates = 117; number of bins = 26; global p value = 0.0037.
ing the Late Pleistocene and Early Holocene in the eastern Adriatic and western Balkans, the SCPD method was used in order to reconstruct population dynamics.The results have confirmed the boom and bust pattern in the Balkan region, and have also shown that these processes happened within the same time frame in the eastern Adriatic (Pilaar Birch, Vander Linden 2017.Figs. 5 and 6).
Additional data are available for the territory of Croatia.Boti≤ presents the results of summed distributions of the Star≠evo (Early Neolithic) and Sopot (Late Neolithic) dates from the territory of Croatia (Boti≤ 2016.17, Fig. 4).A sign of possible population growth similar to the one observed in data from neighbouring regions is present.However, it should be noted that the number of dates from this study is very low (23 dates).
It is interesting to note that the peaks from the three curves (Bulgaria, Hungary, Serbia) coincide (Fig. 4), given the usual assumption that the Neolithic gradually spread from south to north.The earliest Neolithic in Serbia and Bulgaria is dated to ~6200 and ~6100-6050 cal BC, respectively (Whittle et al. 2002;Krauß et al. 2014).The earliest Körös sites in Hungary are not older than 6000 cal BC (Anders, Siklósi 2012.153).Therefore, it should be expected that population boom in Central and Eastern Balkans should have happened earlier than in the Great Hungarian Plain if the demographic process was the same.Given the small samples in all regions of Southeast Europe and the fact that the SCPD method is a very rough tool, these contradictions should not be given too much weight at this moment, as the precision to discriminate between the shifts of one or two centuries may be lacking in this case.
Given the increasing importance of the research focusing on the relationship between climate changes and cultural dynamics (e.g., Wirtz, Lemmen 2003;Budja 2007;2015;Gronenborn 2009;Weninger et al. 2009;2014;Clare, Weninger 2010;Shennan et al. 2013;Lemmen, Wirtz 2012;Boti≤ 2016;Pilaar Birch, Vander Linden 2017), we compared the SCPD curves to a global climate proxy.In order to explore the relationship between climate and population dynamics in the three regions of Southeast Europe, the SCPD curves for Serbia, Hungary and Bulgaria are plotted against the GISP2 core curve (Grootes, Stuiver 1997) (Fig. 4).There is some indication that the troughs after ~5500 cal BC on SCPD curves based on dates from Serbia and Hungary may be related to the reduction in temperature which started somewhat earlier, but no strong sign of covariation is present.Unfortunately, at this point, no high-resolution climate proxies are available for the study region, and a more precise climate reconstruction is not possible.Relying on global proxies such as GISP cores when investigating populations dynamics in the central Balkan area can only be regarded as a general framework for further research.
In order to create a more accurate and precise reconstruction of population dynamics in Southeastern Europe, it is necessary to generate a new sample of radiocarbon dates according to a probabilistic sampling design, which would minimise the re- search bias and increase statistical power.For the region of Central Balkans, the collection of new radiocarbon samples is currently under way by the authors of this paper.This sample is specifically designed for the purposes of population dynamics reconstruction with the SCPD method; an effort is made to approximate a random sample.Therefore, the results presented in this paper should be considered as preliminary and will be further refined when the new data arrive.

Conclusions
It can be concluded that in the Eastern Balkans and Great Hungarian plain, the shape of the SCPD curve is consistent with the population boom predicted by the NDT theory and empirical results from other parts of Europe and Balkans.No clear influence of climate on population dynamics patterns during the 6 th millennium cal BC was detected.
This research is a result of the Project "BIRTH: Births, mothers and babies: prehistoric fertility in the Balkans between 10 000-5000 BC", funded by the European Research Council (ERC) (https://erc.europa.eu/) under the European Union's Horizon 2020 research and innovation programme (Grant Agreement No. 640557;Principal Investigator: SS).We would like to thank Jugoslav Pendi≤ for his help with Figure 1.We are also grateful to Carsten Lemmen and the anonymous reviewer for their constructive comments, suggestions and criticism.The responsibility for all remaining omissions and errors is exclusively ours.5221 Görsdorf, Bojad/iev 1996.127-128