Hostname: page-component-cd9895bd7-gxg78 Total loading time: 0 Render date: 2024-12-23T11:10:43.241Z Has data issue: false hasContentIssue false
  • English
  • Français

SINE QUA NON: INFERRING KODJADERMEN-GUMELNIȚA-KARANOVO VI POPULATION DYNAMICS FROM AGGREGATED PROBABILITY DISTRIBUTIONS OF RADIOCARBON DATES

Published online by Cambridge University Press:  03 March 2023

Gabriel M Popescu Open the ORCID record for Gabriel M Popescu [Opens in a new window] ,
Cristina Covătaru ,
Ionela Opriș ,
Adrian Bălășescu ,
Laurent Carozza ,
Valentin Radu ,
Constantin Haită ,
Tiberiu Sava Open the ORCID record for Tiberiu Sava [Opens in a new window] ,
C Michael Barton  and
Cătălin Lazăr Open the ORCID record for Cătălin Lazăr [Opens in a new window]
Gabriel M Popescu
Affiliation:
Research Institute of the University of Bucharest, Division of ArchaeoSciences, University of Bucharest, No. 90, Panduri Street, Sector 5, 050663, Bucharest, Romania School of Complex Adaptive Systems, Arizona State University, 1031 S. Palm Walk, 85281-2701, Tempe, AZ, USA
Cristina Covătaru
Affiliation:
Research Institute of the University of Bucharest, Division of ArchaeoSciences, University of Bucharest, No. 90, Panduri Street, Sector 5, 050663, Bucharest, Romania
Ionela Opriș
Affiliation:
Research Institute of the University of Bucharest, Division of ArchaeoSciences, University of Bucharest, No. 90, Panduri Street, Sector 5, 050663, Bucharest, Romania
Adrian Bălășescu
Affiliation:
Research Institute of the University of Bucharest, Division of ArchaeoSciences, University of Bucharest, No. 90, Panduri Street, Sector 5, 050663, Bucharest, Romania Institute of Archaeology “Vasile Pârvan”, No. 11, Henri Coandă Street, Sector 1, 010667, Bucharest, Romania
Laurent Carozza
Affiliation:
Research Institute of the University of Bucharest, Division of ArchaeoSciences, University of Bucharest, No. 90, Panduri Street, Sector 5, 050663, Bucharest, Romania CNRS, UMR 5602 Géode, Université Toulouse 2 Jean Jaurès, Toulouse, France
Valentin Radu
Affiliation:
Research Institute of the University of Bucharest, Division of ArchaeoSciences, University of Bucharest, No. 90, Panduri Street, Sector 5, 050663, Bucharest, Romania National Museum of History of Romania, No. 12, Calea Victoriei, Sector 3, 030026, Bucharest, Romania
Constantin Haită
Affiliation:
National Museum of History of Romania, No. 12, Calea Victoriei, Sector 3, 030026, Bucharest, Romania
Tiberiu Sava
Affiliation:
IFIN-HH, 30 Reactorului St., Măgurele, Ilfov, 077125, Romania
C Michael Barton
Affiliation:
School of Human Evolution and Social Change, Arizona State University, 900 S. Cady Mall, 85287-2402, Tempe, AZ, USA School of Complex Adaptive Systems, Arizona State University, 1031 S. Palm Walk, 85281-2701, Tempe, AZ, USA
Cătălin Lazăr*
Affiliation:
Research Institute of the University of Bucharest, Division of ArchaeoSciences, University of Bucharest, No. 90, Panduri Street, Sector 5, 050663, Bucharest, Romania
*
*Corresponding author. Email: [email protected]
Rights & Permissions [Opens in a new window]

Abstract

Past human population dynamics play a key role in integrated models of understanding socio-ecological change over time. However, little analysis on this issue has been carried out for the prehistoric societies in the Lower Danube and Eastern Balkans area. Here, we use summed probability distributions of radiocarbon dates to investigate potential regional and local variation population dynamics. Our study adopts a formal model-testing approach to the fifth millennium BC archaeological radiocarbon record, performing a region-wide, comparative analysis of the demographic trajectories of the area along lower Danube River. We follow the current framework of theoretical models of population growth and perform global and regional significance and spatial permutation tests on the data. Specifically, we investigate whether populations on both sides of the Danube follow a logistic pattern of steady growth, followed by a major decline over time. Finally, our analysis of local-scale growth investigates whether considerable heterogeneity or homogeneity within the region may be observed over the time span considered here. The results show both similarities and differences in the population trends across the area. Our findings are showcased in relation to the cultural characteristics of the region’s 5th millennium BC societies, and future research directions are also suggested.

Keywords

Balkansfifth millennium BCKGK-VI population dynamicsLower DanubeMiddle Holocene
Type
Research Article
Information
Radiocarbon , Volume 65 , Issue 2 , April 2023 , pp. 463 - 484
Check for updates
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2023. Published by Cambridge University Press for the Arizona Board of Regents on behalf of the University of Arizona

INTRODUCTION

The so-called “Kodjadermen-Gumelniţa-Karanovo VI cultural complex” (KGK-VI) is a component of the Southeastern Eneolithic Block (SEB), defined based on material culture by the cultural-historical archaeologists from Bulgaria and Romania. It was identified within a broad region of the Eastern Balkans and Lower Danube Valley, delimited to the north by the Carpathians mountains, but also the steps area from north of the Danube Delta to the east by the Black Sea and reaching as far as the Rhodope mountains and Olt River to the west and the Aegean in the south. The general chronological position of KGK-VI evolution is placed during most of the 5th millennium BC. It is characterized by the widespread development of tell settlements, the emergence of extramural cemeteries (in some cases including wealthy graves), the advent and development of copper and gold metallurgy, along with consistent changes in lithic and ceramic technologies (e.g., graphite and gold pottery painting, exploitation, and processing of various pigments, etc.). Local variants emerged from a common, prior cultural background (e.g., “Varna Culture” on the western coast of the Black Sea) (Todorova and Zhelyaskova Reference Todorova and Zhelyaskova1978; Todorova Reference Todorova1986; Marinescu-Bîlcu Reference Marinescu-Bîlcu2001; Popovici Reference Popovici2010; Petrova Reference Petrova2016; Lazăr et al. Reference Lazăr, Mărgărit and Radu2018; Chapman Reference Chapman2020) and also in areas of interaction with other cultures (e.g., known as Stoicani-Aldeni/Aldeni II cultural group in north-eastern Wallachia) at the intersection between KGK-VI and Cucuteni-Tryplie communities (Frînculeasa Reference Frînculeasa2016), or in previously unoccupied territories (e.g., known as Bolgrad/Bolgrad-Aldeni group) in the region north of the Danube Delta (Todorova and Zhelyaskova Reference Todorova and Zhelyaskova1978; Todorova Reference Todorova1986; Frînculeasa Reference Frînculeasa2016), which marks the most northern and eastern extension of the KGK-VI.

Moreover, the increased mobility of the KGK-VI population is highlighted by a vast exchange network of raw materials and goods. Thus, pottery painted with black, white, and red pigments, graphite or even with gold, along with metal or flint items (e.g., super-blades), adornments made of Mediterranean shells (e.g., Spondylus sp., Dentalium sp., etc.), along with non-local minerals and pigments (e.g., malachite, marble, carnelian, agate, hematite, kaolin, etc.) are common discoveries in the KGK-VI sites (e.g., tells, flat settlements, off-tells, and cemeteries) located in areas that usually lack all of these raw materials (Bailey Reference Bailey2000; Todorova Reference Todorova and Grammenos2003; Anthony et al. Reference Anthony and Chi2010; Popovici Reference Popovici2010). All of this reflects the adoption of new lifeways, economic development and a new way of environmental exploitation and control, in parallel with the development of complex and stratified society and the emergence of elites, as demonstrated by the wealthy graves from Varna I cemetery or other domestic/funerary discoveries (Bailey Reference Bailey2000; Chapman et al. Reference Chapman, Higham, Slavchev, Gaydarska and Honch2007; Anthony et al. Reference Anthony and Chi2010; Slavchev Reference Slavchev2010; Chapman Reference Chapman2020). Consequently, starting from that flourishing development of the SEB human groups, without any parallel in the past, some scholars described the period as the “Ex Balcanae Lux” phenomenon (Todorova and Zhelyaskova Reference Todorova and Zhelyaskova1978; Sterud et al. Reference Sterud, Evans and Rasson1984), the “Climax Copper Age” (Chapman et al. Reference Chapman, Higham, Slavchev, Gaydarska and Honch2007) or the “Golden 5th Millennium BC” (Boyadzhiev and Terzijska-Ignatova Reference Boyadzhiev and St2011), in order to highlight the extraordinary progress of the Neolithic communities.

Within this large region, different human groups have developed specific material culture signatures, or archaeological “cultures,” during this period. Recent aDNA investigations in southeastern Europe have demonstrated that populations responsible for the above mentioned “cultures” of the second half of the 6th millennium and 5th millennium BC, KGK-VI included, have similar genetic features and common ancestry due to their common origin in southwestern Anatolia (Hervella et al. Reference Hervella, Rotea, Izagirre, Constantinescu, Alonso, Ioana, Lazăr, Ridiche, Soficaru and Netea2015). Along with the Anatolian Neolithic ancestry, those populations showcase sporadic evidence of steppe-related ancestry (specimens in Varna I and Smyadovo cemeteries in Bulgaria), but also a consistent hunter-gatherers related ancestry (some resilient native Mesolithic groups from the target area), which indicates a complex population structure and admixture, with several genetic components (Mathieson et al. Reference Mathieson, Alpaslan-Roodenberg, Posth, Szécsényi-Nagy, Rohland, Mallick, Olalde, Broomandkhoshbacht, Candilio and Cheronet2018).

The last decade has seen major advances in developing theoretical, analytical, and methodological instruments, concerning the understanding of demographic change out of large datasets of archaeological radiocarbon (14C) dates, in different parts of the world and encompassing large temporal spans of the Late Pleistocene and Holocene (Shennan et al. Reference Shennan, Downey, Timpson, Edinborough, Colledge, Kerig, Manning and Thomas2013; Downey et al. Reference Downey, Bocaege, Kerig, Edinborough and Shennan2014; Timpson et al. Reference Timpson, Colledge, Crema, Edinborough, Kerig, Manning, Thomas and Shennan2014; Crema et al. Reference Crema, Habu, Kobayashi and Madella2016; Downey et al. Reference Downey, Haas and Shennan2016; Bevan et al. Reference Bevan, Colledge, Fuller, Fyfe, Shennan and Stevens2017; Barton et al. Reference Barton, Aura Tortosa, Garcia-Puchol, Riel-Salvatore, Gauthier, Vadillo Conesa and Pothier Bouchard2018; Riris Reference Riris2018; Roberts et al. Reference Roberts, Woodbridge, Bevan, Palmisano, Shennan and Asouti2018; Timpson et al. Reference Timpson, Barberena, Thomas, Méndez and Manning2020; Crema and Shoda Reference Crema and Shoda2021). Some of these studies have also been focused on understanding the population dispersals and change during the Neolithic of Southeastern Europe, although mostly on the early Neolithic of the region (Porčić and Nikolić Reference Porčićić and Nikolić2016; Blagojević et al. Reference Blagojević, Porčić, Penezić and Stefanović2017; Harper Reference Harper2019; Vrhovnik Reference Vrhovnik2019; Porčić et al. Reference Porčić, Blagojević, Pendić and Stefanović2021; Vander Linden and Silva Reference Vander Linden and Silva2021).

The current study aims at contributing to the understanding of the population geo-temporal dynamics of one of the archaeological signals of the Chalcolithic in Eastern Europe, namely KGK-VI, which reflects the maximum point of development of human communities that lived a Neolithic way of life in this part of Europe. More specifically, our research explores (1) the nature of population trajectories within the chosen region, on both sides of the Danube and (2) the extent to which KGK-VI differs north and south of Danube. The analyses used in the study allow for local and global tests of significance to be performed and regional population histories to be compared through the comparison of empirical and simulated summed probability distributions (SPDs) of radiocarbon dates (see Bevan et al. Reference Bevan, Colledge, Fuller, Fyfe, Shennan and Stevens2017; Crema and Bevan Reference Crema and Bevan2021 for details).

DATA

In the current study 440 radiocarbon dates from both Romania and Bulgaria ascribed to the KGK-VI were used (Figure 1). The data were acquired from publications, gray literature, unpublished sources, and from on-line databases (Reingruber and Thissen Reference Reingruber and Thissen2017; Weninger et al. Reference Weninger, Joris and Danzeglocke2020) (e.g., http://www.14sea.org/index.html and https://www.academia.edu/40774947/CalPal_Holocene_Palaeolithic_14C_Database). The available radiocarbon dates from the above-mentioned sources were compiled in a database (n = 257 from Romania and n = 183 from Bulgaria), grouped into 135 bins for analysis (see explanation below of binning), recovered from 59 sites (n = 32 in Romania and n = 27 in Bulgaria). The database includes contextual information such as site name, site id, site recovery context, “culture” and phase (where available), region, country, laboratory number, the uncalibrated date and uncertainty, and geographical coordinates (longitude and latitude) (Supplementary Table 1).

Figure 1 Map of radiocarbon distribution data set. Legend: 1-Akladi Cheiri; 2-Bikovo; 3-Čardako-Slatino; 4-Djakovo; 5-Dolnoslav; 6-Drama-Merdžumekja; 7-Durankulak; 8-Ezero; 9-Goljamo Delčevo; 10-Hotnica; 11-Junacite; 12-Karnobat; 13-Košarna; 14-Omurtag; 15-Orlitsa; 16-Ovčarovo; 17-Povelyanovo; 18-Smjadovo; 19-Sušina; 20-Tatul; 21-Tell Azmak; 22-Tell Karanovo; 23-Tell Russe; 24-Varhari; 25-Varna1; 26-Varna2; 27-Varna3; 28-Baia Boruz Tell; 29-Popina Blagodeasca; 30-Bordușani; 31-Carcaliu; 32-CăscioareleOstrovel; 33-Cunești; 34-Dambul lui Haralambie; 35-Gumelnița-terrasse; 36-Gumelnița-tell; 37-Hârșova; 38-Lișcoteanca-Movila Olarului; 39-Lunca; 40-Luncavița; 41-Mălăieștii de Jos; 42-Măriuța-C; 43-Măriuța-T; 44-Navodari; 45-Niculițel; 46-Orbeasca Sus; 47-Panduru; 48-Pietrele; 49-Seciu; 50-Șeinoiu; 51-SultanaGhețărie; 52-Sultana-Malu Roșu-terrasse; 53-Sultana-Malu Roșu-tell; 54-Taraschina; 55-Taraschina_2; 56-Urlați; 57-Vărăști; 58-Vitănești; 59-Vlădiceasca.

The dates selected to use in the current investigation were obtained from samples from various materials including, wood, charcoal, seeds, and bones (herbivores and humans). The dates based on shell, fish and carnivore samples have not been included to avoid potential reservoir effects issues. We also applied a cleaning protocol and excluded all dates with large uncertainties in order to remove any potential spurious dates. However, most of our dates have very low uncertainties, with a mean of 43 and a median of 37 years. Recent studies, dedicated to similar research agenda, have also shown that an overly strict data cleaning and the exclusive use of dates with very low uncertainties might potentially be just as damaging for this kind of analysis just as an uncritical data collection (Timpson et al. Reference Timpson, Colledge, Crema, Edinborough, Kerig, Manning, Thomas and Shennan2014, Reference Timpson, Manning and Shennan2015; Vander Linden and Silva Reference Vander Linden and Silva2021).

METHODS

Our analysis was carried out with the R environment for statistical analysis (R Core Team 2020), and the rcarbon R-package for date calibration and SPD modeling (Bevan et al. Reference Bevan, Crema, Bocinsky, Hinz, Riris and Silva2020; Crema and Bevan Reference Crema and Bevan2021), using the Northern Hemisphere calibration curve (IntCal20) (Reimer et al. Reference Reimer, Austin, Bard, Bayliss, Blackwell, Ramsey, Butzin, Cheng, Edwards and Friedrich2020), along with gstat (Pebesma and Graeler Reference Pebesma and Graeler2021), sp (Pebesma et al. Reference Pebesma, Bivand, Rowlingson, Gomez-Rubio, Hijmans, Sumner, MacQueen, Lemon, Lindgren and O’Brien2022), and other R packages mentioned in the rmarkdown script that accompanies the study. Dataset, supplemental figures and R Markdown scripts needed for reproducing the results of our analysis are available at: https://zenodo.org/record/7587242#.Y9hI9S8Ro4c. The methods that we used in our analysis are based on previously developed quantitative analysis of SPDs by Shennan et al. (Reference Shennan, Downey, Timpson, Edinborough, Colledge, Kerig, Manning and Thomas2013), further refined in several other recent studies (Downey et al. Reference Downey, Haas and Shennan2016; Brown and Crema Reference Brown and Crema2019; Crema et al. Reference Crema, Habu, Kobayashi and Madella2016; Riris Reference Riris2018; Timpson et al. Reference Timpson, Colledge, Crema, Edinborough, Kerig, Manning, Thomas and Shennan2014; Timpson et al. Reference Timpson, Manning and Shennan2015), the non-parametric extension devised by Crema et al. (Reference Crema, Habu, Kobayashi and Madella2016; see also Crema and Kobayashi Reference Crema and Kobayashi2020), and the spatial permutation test (Crema et al. Reference Crema, Bevan and Shennan2017; Crema and Bevan Reference Crema and Bevan2021). The reader should refer to these bibliographical resources (and references therein) for more details on the methods and the concepts behind them.

Spatial Dispersal Analysis

We first analyzed the spatial structure of the dispersal of our 14C distribution (Figure 2a,b). We considered this as being the first step in our analysis for a better understanding of the spatial dynamics of the dispersal of the KGK-VI to assess its correlation (potentially) with demography. Based on Hengl (Reference Hengl2006) equation for establishing the right pixel size for a grid, and our regional context we superimposed a 23 × 23 km grid on the area, with each grid cell covering approximately 460 km2. The analysis was done in R statistical environment with the gstat (Pebesma and Graeler Reference Pebesma and Graeler2021) and sp (Pebesma et al. Reference Pebesma, Bivand, Rowlingson, Gomez-Rubio, Hijmans, Sumner, MacQueen, Lemon, Lindgren and O’Brien2022) packages. We divided the study area into grid cells and hexagon cell shapes were chosen, given their shape being the closest to a circle and the easy to use in a tessellation. Their minimal edge effects as well as the identical neighboring cells and having the same distance between centers for all the neighbors, make them particularly suitable for our analysis (see also Vrhovnik Reference Vrhovnik2019).

Figure 2 (a) Number of radiocarbon dates per grid cell. Values are log10 scaled. (b) The earliest appearance of the KGK-VI settlements in the grid cells, consisting of grid cell centroids with the date for the beginning of the KGK-VI occupation. Gridded area (white hexagons) represents the KGK-VI area with dated sites.

Thus, the study area is covered with 477 grid cells, of which only 47 grid cells are occupied with sites, forming several clusters, and a patchy distribution of samples. The number of radiocarbon dates per grid cell varies from one (seen in 11 grid cells) to a maximum of 67 (seen in one cell) with a median value of 3 dates per grid cell and third quartile at 10.5 dates per grid cell. The distribution of dates per cell is shown in Figure 2a.

For each grid cell, we then calculated a normalized summed calibrated radiocarbon probability distribution. To calculate the calendar age ranges, highest probability density was used, and these are the shortest ranges that include 95% of the probability in the summed probability density function. As such, the starting date of the KGK-VI in a particular grid cell was taken to be the lower 95% range endpoint date. These estimated starting dates are shown in Figure 2b. Consequently, these dates were to estimate the spread of the KGK-VI across the area. Grid cells with only one radiocarbon date were excluded from the interpolation.

Dates Binning, Calibration, and SPDs Production

To mitigate potential issues due to the differences in intensity of sampling, radiocarbon dates are combined in 100-year bins within each archaeological context (e.g., horizontal and vertical provenience units) so that the intensively sampled sites/areas are not overrepresented and cause artificial spikes in the observed SPDs. When multiple dates are present from a single site, they are aggregated within each archaeological context and when their distance in 14C years is less than 100 years. The procedure applies a hierarchical cluster analysis using the complete linkage method and a cut-off value of 100 years to separate the observations. Although our selection of 100 years for binning is arbitrary, we performed a bin sensitivity analysis (Supplementary Figure 1), which shows this choice has no negative impact on the accuracy of results (all bin sizes fit within the 95% confidence simulated envelope—gray area) and is also well above the median error of 37 years in the dataset (our protocol follows those already applied in the literature; for more details (Bevan et al. Reference Bevan, Colledge, Fuller, Fyfe, Shennan and Stevens2017; Riris Reference Riris2018; Crema and Bevan Reference Crema and Bevan2021).

Dates were calibrated using the Northern Hemisphere Radiocarbon Age Calibration Curve (IntCal20) (Reimer et al. Reference Reimer, Austin, Bard, Bayliss, Blackwell, Ramsey, Butzin, Cheng, Edwards and Friedrich2020) and the rcarbon package (Bevan et al. Reference Bevan, Crema, Bocinsky, Hinz, Riris and Silva2020). Multiple dates within a bin are calibrated and summed “inside” the bin and subsequently divided by the number of dates so that each archaeological context contributes a single date distribution to the overall SPD. The probability distributions of the calibrated dates were summed over the entire KGK-VI period to produce empirically based SPDs using the entire data set, as well as for subsets for the two regions north and south of the Danube. Following detailed discussions in recent works (Weninger et al. Reference Weninger, Clare, Jöris, Jung and Edinborough2015; see also details in Bevan et al. Reference Bevan, Colledge, Fuller, Fyfe, Shennan and Stevens2017), we have not normalized the post-calibration distribution of each date (that ensures it sums to 1 under the curve) before summation of multiple dates. This ensures the reduction of creation of abrupt spikes in the final summed probability distributions, there where the calibration curve is steep (Weninger et al. Reference Weninger, Clare, Jöris, Jung and Edinborough2015; Bevan et al. Reference Bevan, Colledge, Fuller, Fyfe, Shennan and Stevens2017; Crema and Bevan Reference Crema and Bevan2021).

Model Testing

In order to differentiate SPD fluctuations that represent meaningful demographic change from those due to sampling error noise, we compare the observed SPDs against null models of simulated dates derived from hypothesized calendar age distributions of dates. These null (hypothesized) models assume increasing survival of datable radiocarbon material through time according to either an exponential or a logistic population growth. Only portions of the SPD curves that fall outside the 95% confidence interval (CI) of the null models are considered sufficiently significant demographic changes to be considered in our subsequent discussion. Model testing procedures are described in detail below.

We first evaluate the goodness-of-fit of the entire data set SPD, by first fitting the calibrated data to a generalized exponential model, with the help of modelTest function (Figure 3) (Timpson et al. Reference Timpson, Colledge, Crema, Edinborough, Kerig, Manning, Thomas and Shennan2014; Crema and Bevan Reference Crema and Bevan2021) (Supplemental Material for analysis in R). An exponential model had become common practice, as it is assumed to account for population growth with unlimited resources and taphonomic processes. We assessed whether the SPDs of the 14C dates for the entire region showed statistically relevant deviations when compared against the exponential model, following the procedure described in several other studies (Timpson et al. Reference Timpson, Colledge, Crema, Edinborough, Kerig, Manning, Thomas and Shennan2014; Bevan et al. Reference Bevan, Colledge, Fuller, Fyfe, Shennan and Stevens2017; Crema and Bevan Reference Crema and Bevan2021; Crema et al. Reference Crema, Habu, Kobayashi and Madella2016; Riris Reference Riris2018; Palmisano et al. Reference Palmisano, Bevan and Shennan2017; Roberts et al. Reference Roberts, Woodbridge, Bevan, Palmisano, Shennan and Asouti2018 and references therein).

Figure 3 Results of fitting and comparing the entire regional empirical SPD against the exponential null model of population growth. Monte Carlo 95% confidence null model gray envelope is based on 5000 runs. Observed SPD is shown with solid red line, while the positive and negative deviations from the null are marked in red and blue. (Please see online version for color figures.)

Null models were simulated for the entire KGK-VI region using assumptions of the exponential and logistic growth patterns in the following way, repeated for each type of model.

  1. 1. An exponential/logistic model was generated for the entire time period and fit to the empirical SPD produced by the dates we compiled.

  2. 2. A set of dates (equivalent in number to the bins used for the empirical SPD) was generated from the model and errors assigned randomly (within the range of empirical date errors).

  3. 3. An SPD was generated from the model dates and errors.

  4. 4. Steps 2–3 were repeated 5000 times to estimate the 95% CI around the model.

The empirical SPDs are then compared with the 95% CI around each of the two models (exponential and logistic). Portions of the observed regional SPD that fall outside the envelope were considered statistically significant local deviations above and below the null model (red and blue areas respectively). Following methods outlined by Timpson et al. (Reference Timpson, Colledge, Crema, Edinborough, Kerig, Manning, Thomas and Shennan2014) a global p-value can then be calculated from the total area of the empirical SPD curve that falls outside the 95% CI of each null model.

To evaluate how well the exponential model fit to our empirical data we employed Akaike Information Criterion (AIC) to select the most parsimonious fitted model for the entire region (Sakamoto et al. Reference Sakamoto, Ishiguro and Kitagawa1986). AIC suggested that the use of different model might be a better fit for our context. As such, we decided to use the logistic growth model as the null model for comparison and discussion of results (Figure 4), following the same procedures outlined above for the exponential growth model. We also used an extension of the global p-value to evaluate the point-to-point differences along the SPD curve. This test compares the empirically observed difference to the distribution of differences in the SPD curve between two points in time against the distribution of expected values under the null hypothesis (Edinborough et al. Reference Edinborough, Porčić, Martindale, Brown, Supernant and Ames2017). We also deployed a more recent alternative to SPDs, Composite Kernel Density Estimates (CKDEs) (Brown Reference Brown2017; McLaughlin Reference McLaughlin2019), which has the advantage of minimizing calibration artificial spikes, as well as, providing estimates of sampling and calibration-derived uncertainty over time (Figure 5). Based on the qualitative inspection of the SPDs, the CKDE curve and formal AIC test, we decided to also fit a suite of four composite models (see Rmarkdown document in the Zenodo repository) to the summed calibrated probability distribution for KGK-VI, representing potential demographic models (see Goldberg et al. Reference Goldberg, Mychajliw and Hadly2016; Arroyo-Kalin and Riris Reference Arroyo-Kalin and Riris2021; de Souza and Riris Reference de Souza and Riris2021). Assessing for the goodness of fit of these models was achieved following the same procedure as outlined above.

Figure 4 Results of fitting and comparing the entire regional empirical SPD against the logistic null model of population growth. Monte Carlo 95% confidence null model gray envelope is based on 5000 runs. Observed SPD is shown with solid red line, while the positive and negative deviations from the null are marked in red and blue.

Figure 5 Bootstrapped composite kernel density estimate, suggesting a composite model with the breakpoint at approximately 4400 BC should be tested.

Permutation Tests—Regional Comparison

We are obviously interested in empirically testing the variation between the northern and southern regions of the KGK-VI (Danube River is used as a geographical divide). In order to achieve this, each region was compared to a null model assuming no spatial differences. This null model can be obtained by pooling the radiocarbon dates from across the entire region and simulating from the pooled SPD (Figure 6).

Figure 6 The KGK-VI empirical SPD record fitted to logistic (5000 BC–4400 BC) and exponential (4400 BC–5750 BC) models, with significance envelope derived from 5000 Monte Carlo simulations.

Crema et al. (Reference Crema, Habu, Kobayashi and Madella2016) developed a permutation-based test to statistically compare two or more SPDs (Figure 7). The null hypothesis is generated by simulating multiple (e.g., 5000 here) SPDs whose dates are drawn randomly from both subregions (north and south of the Danube). These simulated SPDs are again combined, as above, to produce a 95% CI envelope against which the SPD from each region can be compared (see Crema et al. Reference Crema, Habu, Kobayashi and Madella2016; Crema and Bevan Reference Crema and Bevan2021).

Figure 7 Permutation test showing variation between regional population growth. Observed SPDs for each region are shown with a solid black line, while the dashed line represents the observed pan-regional SPD. Gray areas represent the 95% confidence envelope for the null model, red and blue bands represent areas where the observed SPD significantly positively (red) and negatively (blue) deviates from the pan regional null model.

This is a robust approach to inter-regional differences in the research intensity because the comparison is based on the shape of the SPD (the relative change in summed probabilities within each region) and not on differences in their absolute magnitudes. Maintaining the observed number of bins for each region and comparing population trajectories rather than absolute differences in density, the permutation test bypasses the problem. Moreover, sample size is taken into account in the width of the 95% CI envelope. Significant negative (or positive) deviations of the SPD in one region does not necessarily imply a lower (or higher) absolute population density, but that the drop in the proxy within the dynamics of that region was significantly stronger compared to rest of the data.

Spatial Permutation Test

The spatial permutation test is an extension of the permutation test described above, having the virtue of allowing for the assessment of variation without the imposition of a priori regions of analysis. The steps involved in the spatial permutation test are described in detail by Crema et al. (Reference Crema, Bevan and Shennan2017; see also Crema and Bevan Reference Crema and Bevan2021). Below are summarized the steps involved in the spatial permutation analytical protocol.

  1. 1. Produce local SPDs for each site with dates combining the date distribution at the site with date distributions at neighboring sites weighted as a function of their distance from the site. We selected a neighborhood radius of 100 km following a sensitivity analysis of different radii (see Supplemental Figures 2–3).

  2. 2. Divide the KGK-VI temporal span of 5050–3800 BC into equal transition blocks (e.g., time slices), 250 years each, which, given the length of our time span, we consider relevant in our endeavor to detect long term regional growth patterns (Figure 8).

    Figure 8 Observed rate of growth at each transition computed from the SPD. I: 5050–4800 to 4800–4550; II: 4800–4550 to 4550–4300; III: 4550–4300 to 4300–4050; IV: 4300–4050 to 4050–3800 BC.

  3. 3. Calculate the overall growth rate (change in the SPD curve) between each temporal transition block and the subsequent one.

This allowed us to evaluate spatial patterns of demographic growth and decline in two ways. We compared the growth rates calculated for each local SPD with the overall growth rate for each transition block (Figure 9). For each transition block, we also compared the growth rates of local SPDs at each site with a simulated model generated by repeatedly (10,000 iterations) randomly shuffling the local SPDs spatially across all site locations and combining the growth rates of the shuffled SPDs at each site. This allowed us to identify “hot” and “cold” spots (areas of significance), defined as areas where the local growth exceeds the growth observed in the simulation (Figure 10). Following methods discussed in Crema et al. (Reference Crema, Bevan and Shennan2017), two measures of significance are produced in the course of the spatial permutation test. p-values are measures of significance between observed local growth and simulated growth rates. However, the use of multiple testing approach, increases the potential for compounding false positive results (e.g., some local SPDs will be higher or lower than the theoretical expectation by chance alone). A more robust q-values test is therefore also computed by adjusting p-values to account for false positive discovery rate. Thus, a p-value of 0.05, implies that 5% of the tests will result in false positives, a q-value of 0.05 means that 5% of the results that have a q-values less than 0.05 are false positives (see Crema et al. Reference Crema, Bevan and Shennan2017 for further details).

Figure 9 Local geometric growth rate for each transition block. I: 5050–4800 to 4800–4550; II: 4800–4550 to 4550–4300; III: 4550–4300 to 4300–4050; IV: 4300–4050 to 4050–3800 BC; V: Geographical reference map for the local geometric growth rate analysis, shown in transitional blocks I-IV.

Figure 10 Spatial permutation test showing where growth is significantly higher or lower than the null for each transition block. I:5050–4800 to 4800–4550; II:4800–4550 to 4550–4300; III:4550–4300 to 4300–4050; IV:4300–4050 to 4050–3800 BC; V: Geographical reference map for the spatial permutation test analysis, shown in transitional blocks I-IV. Significance is shown in terms of q-values (more robust against false positives) and p-values.

RESULTS

Spatial Dispersal Analysis

Figure 2(a–b) displays the results of spatial interpolation based on the analysis of the earliest appearance of the KGK-VI settlements in the grid cells. Early presence of KGK-VI across the region, between ∼4900 BC–4700 BC, although rather patchy, seems to be documented in several spots across the region, but mostly in the central and southern KGK-VI area, along major river valleys and their surroundings (see also in conjunction with Figure 1). The interpolation results thus suggest a fast emergence of this archaeological signature on both sides of Danube River, in its main core areas, in central-northeastern Bulgaria and southern Romania, extending in southeastern Romania, southern and southwestern Bulgaria (Figure 2b). There also seems to be a lag in KGK dispersal, roughly around 4600–4550 BC, or a weaker signal, that could be also due to a smaller sample of dated sites. From these core area KGK-VI dispersed in most of the region by 4400 BC. This analysis also suggests that potentially two more clusters, seems to emerge, one in south-southwestern Bulgaria and the other in southeastern Romania. Given the number of sites with radiocarbon dates, available for our study, we have not split the database into more clusters, as this might be the result of the current sample size and might have created artificial patterns in the SPD analysis and its model testing.

Model Testing

Figure 3 presents three noticeable fluctuations in the empirical SPD, from the exponential null model. The first of those is a significant negative departure from the exponential null model between 4895 BC–4661 BC (marked in blue). This is followed by a long significant positive deviation from the null between 4560 BC and 4281 BC (marked in red), with its peak at ∼4408 BC. A point-to-point statistical test for changes in the shape of the SPD, suggests that the difference between the peak in the empirical SPD and a first major trough at 4175 BC are statistically significant at a p = 0.0004 level. After this highest point the SPD drops rapidly, while remaining above the 95% CI null envelope, until ∼4300 BC. Subsequently, the SPD remains within the null model 95% CI envelope until dropping below this limit after 4200 BC until the end of the sequence, thus marking the second statistically significant negative deviation from the null model. The data, overall, exhibit significant deviations from the exponential model at a global level of p = 0.0002.

The result provided by the AIC test, comparing the fitted exponential and logistic null models indicated the logistic null model as more parsimonious (ΔAIC = 235, 42).

We reran the modelTest function with the logistic null model, and its result are shown in Figure 4. The SPD deviates statistically from the logistic model as well, with an overall p = 0.0002. A short positive deviation from the null model at the beginning of the SPD is most likely a statistical artifact of the method, which fits the model at a later start date. Inspecting Figure 4 we observe a better agreement between the fitted model (dashed line) and the 95% CI envelope, although toward the end, the envelope departs from the fitted line, and is more in line with the pattern of the empirical SPD. At the same time, the empirical SPD is in good agreement with both the null model envelope and the fitted model for about a half of the time range. The beginning and early stages of KGK-VI follow a steady growth pattern within the limits of the null model, until a rapid rise that starts approximately around 4600 BC and significantly deviates positively from the null and fitted line model, from 4553 BC until 4303 BC, reaching the peak around 4400 BC. The result of the point-to-point test that we have run between the highest peak in the SPD and the first major subsequent trough in the SPD around 4180 BC, also proved to be statistically significant at p = 0.0004 level. A final major deviation from the logistic model is represented by a long significant departure from the null, marking the final stages of the KGK-VI.

The bootstrapped CKDE analysis (Figure 5), as well as the shape of the empirical SPDs, and the significantly departure from both the exponential and logistic growth models, confirmed the inadequacy of both models to explain the KGK-VI data. We found that the best fit model (a Logistic-Exponential model (Figure 6), adequately describes the observed SPD by combining a phase of logistic growth (∼5000–4400 cal BC), followed by an almost symmetrical, this time exponential decrease (∼4400–3800/3750 cal BC), with a global p-value of 0.5257 (see also R code in Supplemental Material for the analytical protocol). It is worth mentioning that the empirical SPD, highly matches the 95% confidence envelope of the composite null model and that, while neither local nor global deviation from the null emerges, the decreasing trend after 4400/4350 BC is continuous with no signs of positive peaking, until the end of the KGK-VI archaeological signature.

Permutation Test—Regional Comparison

The permutation approach compares the regional trajectories from the two major areas of the KGK-VI, namely north and south of Danube, with the model representing the general growth trends of the targeted region, as opposed to an idealized model. Figure 7 shows that both regions significantly deviate from the whole region null model at p = 0.001 for the southern area, and p = 0.001 for the northern area, respectively. The overall results, thus suggests regional differences relative to the timing and nature of population change.

The south of the Danube area (Figure 7 bottom) starts positively beyond the limits of the regional null model, and the pan-regional SPD (dashed line) as well, until around 4580 BC, growing faster positively deviating from the null. The growth continues afterward but within the limits of the null model and reaches its peak at ∼4500 BC, but faster than the pan-regional empirical SPD (e.g., south of the Danube demography follows a faster pattern than the entire region as a whole at this time). Demographic stability is maintained at the highest density between approximately 4500–4400 BC, followed by an abrupt declining until 4200 BC, and continuing at a steadier pace until the end of the sequence. Two significant negative deviations from the null occur in the southern SPD between 4359–4320 BC and 4243–4222 BC. The north of the Danube area (Figure 7 top), on the other hand, displays a quite different trend initially displaying a significant lag against the null model, until approximately 4550 BC, followed by a faster growth, but closely matching both the regional null model (gray shaded area), and the pan-regional empirical SPD. Unlike the southern area, the northern density peak is reached at the same time as the region wide SPD (dashed line). The maximum density is maintained for about 100 years, like in the southern area. The declining pattern in the northern SPD follows a similar path like the southern one. The difference here is that there are two significant positive deviations from the null between 4358–4320 BC and 4256–4204 BC, respectively, contemporary with the negative ones in the southern region. Most of the local SPD declining pattern, however, falls within the 95% CI of the pan-regional growth model.

Spatial Permutation Test

The final step in our analysis evaluates the spatial variation in population growth without an a priori regional classification but based on the geographical coordinates (see details in previous section). Figure 8 shows the pattern of the observed pan-regional geometric growth rate and in the four intervals across the five 250-year transition blocks at which growth was measured (I: 5050–4800 BC to 4800–4550 BC; II: 4801–4550 BC to 4550–4300 BC; III: 4550–4300 BC to 4300–4050 BC; IV: 4300–4050 BC to 4050–3800 BC). This shows that the pan-regional trend had an initial rapid growth followed by a declining reduction of growth that extends even below zero once the carrying capacity is reached.

Evaluation of the local growth rates across the four temporal blocks (see Figure 9), shows some regional variation in actual growth rates, especially between the transitional blocks I-II and III-IV. To test the significance of actual growth rates and overcome potential issues caused by the calibration, we compared the local growth patterns to a null model based on regional wide growth trends (Figure 10). This allowed for the identification of significant positive or negative deviations from the overall trend. For the first interval of 5050–4800 BC to 4800–4550 BC (Figure 9), all KGK-VI sites from both sides of the Danube River present positive growth rates. However, as shown in Figure 10, most of the sites do not significantly positively deviate from the overall demographic trend, and only a few of them displaying negative deviation at p < 0.05. In the second temporal interval of 4800–4550 BC to 4500–4300 BC (Figure 9), while mostly displaying positive growth, some regional differences between north and south of Danube are apparent. Growth rates in the south range from 0.0004–0.004% annual rate, as opposed to 0.004–0.009% annual rate in northern area. The significance test (Figure 10) also shows obvious differences between the northern and southern KGK-VI areas. As such, the actual positive growth rates significantly depart from the regional trend in most of the northern region (q < 0.05), while the south trends in the opposite direction, with most of the area showing negative deviation from the overall pattern model, p < 0.05 level (see more details in e.g., Crema et al. Reference Crema, Bevan and Shennan2017).

The third and fourth temporal intervals, 4550–4300 BC to 4300–4050, and 4300–4050 BC to 4050–3800 BC, respectively, show similar negative rates throughout the KGK-VI distribution area (Figure 9). However, as shown in Figure 10, these rates neither positively nor negatively deviate from the expectation of the pan-regional model, across the entire area.

DISCUSSION

We have developed in this research a set of radiocarbon-based demographic models of a well-known Chalcolithic archaeological signature in the Eastern Balkans and Lower Danube area, namely KGK-VI. This work contributes to other recent studies assessing Neolithic demographic dynamics in Southeastern Europe (Porčić et al. Reference Porčić, Blagojević and Stefanović2016; Blagojević et al. Reference Blagojević, Porčić, Penezić and Stefanović2017; Vrhovnik Reference Vrhovnik2019; Vander Linden and Silva Reference Vander Linden and Silva2021).

Using the spatial interpolation of the radiocarbon record we were able to show the spatiotemporal structure of the spread of the KGK-VI archaeological signature. It was marked by an early patchy manifestation in the main core area, as well as a few other regions both north and south of Danube. Thus, the analysis suggests a faster emergence of the KGK-VI archaeological signal in central-eastern, north-eastern areas of current Bulgaria, and spread across the entire studied area Although there is a slight lag in the expansion between south and north, the KGK-VI extends its entire occupation approximately simultaneously by about 4500 BC. However, given it rather rapid emergence and spread in some spots of the KGK-VI area, it is difficult, at this time, and it is beyond of this study aims, to determine whether this was due to a population dispersal, or to a rather demic-cultural combined model. This is something to be further modeled in other studies dedicated to this subject.

The positive rejection of the exponential model between 4560 BC and 4281 BC, demonstrates a higher-than-expected growth during this segment. The negative deviation from a null model in the early phase agrees with such a logistic growth trend. It may also be influenced by the slower pace in the northern area. The declining trend is also at a steeper and faster rate than expected, possibly indicating that population overshot regional carrying capacity. This agrees with demographic trends noted for the Neolithic observed in other areas of Europe (see below). The results of the regional permutation test display the regional differences between the north and south of Danube densities, however mostly in the details of their respective SPDs. As such, the northern area does show a lag in the local SPD, compared to the pan-regional null, while the southern area, displays a faster than expected pace toward reaching the peak in its density. However, looking at the shape of the curve of each region one can distinguish that their overall shape is quite similar and that their peak in summed densities fall within the 95% confidence pan-regional envelope, differing only in their position above (in the south) and partially below (in the north) the pan-regional SPD. The fact that the two regions show reciprocal patterns of demographic expansion and decline, suggest that we are potentially seeing some population shifting from the north to the south side of the Danube ca. 4600/4550 BC, and then some movement back the other way during the period of demographic decline after 4400/4350 BC.

Several studies of Early and Middle Holocene societies in other European regions, have shown “booms” followed by important “busts” (drops) in the SPDs (Shennan et al. Reference Shennan, Downey, Timpson, Edinborough, Colledge, Kerig, Manning and Thomas2013; Bernabeu Aubán et al. Reference Bernabeu Aubán, García Puchol, Barton, McClure and Pardo Gordó2016; Downey et al. Reference Downey, Haas and Shennan2016; Bevan et al. Reference Bevan, Colledge, Fuller, Fyfe, Shennan and Stevens2017; Roberts et al. Reference Roberts, Woodbridge, Bevan, Palmisano, Shennan and Asouti2018; Palmisano et al. Reference Palmisano, Woodbridge, Roberts, Bevan, Fyfe, Shennan, Cheddadi, Greenberg, Kaniewski and Langgut2019), which are proposed to be related to long-term decline of Neolithic farming society. A similar pattern is shown in our study, when after a statistically significant “boom” and an equilibrium (of about 100 years—ca. 4450–4350 BC), a long period of higher than expected decline follows (about 350 years—ca. 4150–3800 BC), statistically significant against both exponential and logistic null models. However, we need a larger database to analyze SPDs for farming and climate and/or environmental co-variates, to help evaluate such potential explanations. Other potential explanations could be related to non-demographic aspects, such as changes in settlement patterns and in settlement size. Small more dispersed villages, during the early stages of KGK-VI, as well as a change from the larger complex tells, off-tell habitations, and flat settlement system from its pinnacle, toward much smaller more dispersed villages, may have increased the number of residential features that were dated, thus inflating the SPDs, and giving the appearance of higher population density, although it did not actually change. This is a research direction that deserves further investigation (see Porčić et al. Reference Porčić, Blagojević, Pendić and Stefanović2021; Supplemental Material S1), as well as using complementary data, such as mortality juvenility index, and its correlation with the shape of the SPDs (see details in Downey et al. Reference Downey, Bocaege, Kerig, Edinborough and Shennan2014).

CONCLUSION

In this exploratory study, we provided a possible reconstruction of the population dynamics of one of the well-known archaeological signals of the Chalcolithic in the Balkans and Lower Danube area, namely KGK-VI.

The current approach highlights the following conclusions:

  1. 1. In the targeted geographical area, the human population density steadily grew for the first approximately 300 years following a logistic growth model until it reached the carrying capacity.

  2. 2. The current study clearly shows the existence of at least two initial cores (north and south of the Danube) where what we call KGK-VI originally appeared in the early 5th millennium BC (4901–4751 cal BC) (Figure 2).

  3. 3. The initial dispersion of the KGK-VI archaeological signal from each of the two cores, headed, for the first ca. 300 years, toward south/southwest on either side of the Danube, to disperse to the rest of the region by 4500–4400 BC, reaching the maximum development, documented until now for KGK-VI.

  4. 4. The growth episode was followed by an equally fast start of the decline, beginning at ca. 4350 BC, without a significant recovery, until it vanished from the archaeological record toward the end of the 5th millennium BC–early 4th millennium BC.

  5. 5. This decline was a continuous process, and it covers approx. 550 years (ca. 4350–3800 BC). The underlying causes behind both characteristics, especially relative to the declining trend, require more in-depth research. Therefore, the challenge for the future is to explore whether the “boom” and “bust” is related to deep-time historical contingencies under multiple factors.

  6. 6. Our approach demonstrated the displacement vectors of KGK VI populations in the approached geographical area in correlation with the distinct temporal moments that marked them.

  7. 7. Last, but not least, the present study demonstrates, once again (if necessary!), the fragility of the notion of archaeological “culture” (or cultural complex, group, aspect or facies) concept (-s), and especially that this Kossinian notion creates a major confusion between the notions of population, civilization and ceramic styles “fashion” (to which unfortunately we are still tributary). The entanglement created by these artificial notions in Balkan archaeology (and not only) led to excessive fragmentation of the prehistoric archaeological landscape. It prevented proper understanding (sometimes) of objective eventful realities of the past.

The current study brings some critical corrections related to the absolute chronology of the KGK-VI communities (humans, not pottery style!), especially regarding the end of this civilization, which cannot be placed around 4200 BC as recently postulated by various scholars, but much later. Why some communities in certain settlements are more resilient than others and the SEB collapse causes remain an open question, which we intend to answer in the future.

Furthermore, the research approach outlined here, must involve in future studies on this topic the identification of potential biases, integration of additional dates, and the correlation with other proxies (e.g., mortality juvenility index, data relative to food production, site and house counts, particular pottery style shift, innovations spreading, climate changes, and paleoenvironment alteration) or other constraints (e.g., possible reservoir effect on 14C dating on humans and other omnivorous species), which are critical for further developments of SPD analysis of ancient demography in the Balkans and Southeastern Europe in order to increase the accuracy of the data included in the analysis. Additionally, future studies should be extended to incorporate all Chalcolithic signals in the region, thus allowing the correlation of their results to those provided by other studies elsewhere in Europe and the Near East. In this way, it will be possible to better understand the bigger picture of the past demography, human behaviors determinants (necessities, choices, availability, and constraints), and specific social and economic mechanisms that may connect human groups at larger spatial scales in the targeted area.

ACKNOWLEDGMENTS

Many thanks to our colleagues Radian Andreescu (deceased), Pavel Mirea, Valentin Parnic, Andrei Soficaru, Done Șerbănescu, and Valentina Voinea for providing samples for AMS radiocarbon dating from the archaeological sites where they excavated. Also, we would like to thank Ovidiu Frujină for his involvement in the sampling process and verifying the references. Special thanks go to Oana Gâză, Maria Ilie, Doru Paceșilă, Gabriela Sava, Corina Simion, and Iuliana Stanciu (IFIN-HH, Romania) for their participation in sample preparation and analyzing of radiocarbon samples. We also want to thank two anonymous reviewers for their helpful comments and suggestions that improved the quality of the paper.

This work was funded by UEFISCDI Romania, through the research grant 351PED/2020 (CALIB-RO), project code: PN-III-P2-2.1-PED-2019–4171.

Last, but not least, we dedicate the current study to the memory of Professor Dragomir Popovici (1952–2019), who made his total contribution to the better knowledge of the KGK-VI civilization through the archaeological excavations and multi-disciplinary research from the Hârşova and Borduşani tell settlements (Romania), including a consistent set of radiocarbon data. Sit tibi terra levis.

SUPPLEMENTARY MATERIAL

To view supplementary material for this article, please visit https://doi.org/10.1017/RDC.2023.6

DATA AND MATERIALS AVAILABILITY

All data needed, supplemental figures and R code for reproducing and evaluating the results provided in this paper are published online at: https://zenodo.org/record/7587242#.Y9hI9S8Ro4c

Additional data related to this paper may be requested from the authors.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

References

REFERENCES

Anthony, DW, Chi, J, New York University, editors. 2010. The lost world of old Europe: the Danube valley; 5000–3500 BC. New York: Institute for the Study of the Ancient World.Google Scholar
Arroyo-Kalin, M, Riris, P. 2021. Did pre-Columbian populations of the Amazonian biome reach carrying capacity during the Late Holocene? Philosophical Transactions of the Royal Society B: Biological Sciences 376(1816):20190715. doi: 10.1098/rstb.2019.0715. Available at https://royalsocietypublishing.org/doi/full/10.1098/rstb.2019.0715.CrossRefGoogle Scholar
Bailey, DW. 2000. Balkan prehistory: exclusion, incorporation and identity. London; New York: Routledge.Google Scholar
Barton, CM, Aura Tortosa, JE, Garcia-Puchol, O, Riel-Salvatore, JG, Gauthier, N, Vadillo Conesa, M, Pothier Bouchard, G. 2018. Risk and resilience in the late glacial: a case study from the western Mediterranean. Quaternary Science Reviews 184:68–84. doi: 10.1016/j.quascirev.2017.09.015. [accessed 2019 Nov 4]. Available at http://www.sciencedirect.com/science/article/pii/S0277379117302615.CrossRefGoogle Scholar
Bernabeu Aubán, J, García Puchol, O, Barton, M, McClure, S, Pardo Gordó, S. 2016. Radiocarbon dates, climatic events, and social dynamics during the Early Neolithic in Mediterranean Iberia. Quaternary International 403:201–210. doi: 10.1016/j.quaint.2015.09.020. [accessed 2019 Oct 21]. Available at http://www.sciencedirect.com/science/article/pii/S1040618215008885.CrossRefGoogle Scholar
Bevan, A, Colledge, S, Fuller, D, Fyfe, R, Shennan, S, Stevens, C. 2017. Holocene fluctuations in human population demonstrate repeated links to food production and climate. Proceedings of National Academy of Sciences USA 114(49):E10524–E10531. doi: 10.1073/pnas.1709190114. [accessed 2019 Nov 21]. Available at http://www.pnas.org/lookup/doi/10.1073/pnas.1709190114.CrossRefGoogle ScholarPubMed
Bevan, A, Crema, E, Bocinsky, RK, Hinz, M, Riris, P, Silva, F. 2020. rcarbon: calibration and analysis of radiocarbon dates. Available at https://CRAN.R-project.org/package=rcarbon.Google Scholar
Blagojević, T, Porčić, M, Penezić, K, Stefanović, S. 2017. Early Neolithic population dynamics in the Eastern Balkans and the Great Hungarian Plain. Documenta Praehistorica 44:18–33.CrossRefGoogle Scholar
Boyadzhiev, Y, St, Terzijska-Ignatova, editors. 2011. The Golden Fifth Millennium: Thrace and its neighbour areas in the Chalcolithic. Proceedings of the international symposium “The Golden Fifth Millennium: Thrace and its neighbour areas in the Chalcolithic”, held in 2009 in occasion of the 70th anniversary of the archaeological excavations of Yunatsite Tell. Pazardzhik, 27–30 October 2009. Sofia: Bard.Google Scholar
Brown, AA, Crema, ER. 2019. Māori population growth in Pre-contact New Zealand: regional population dynamics inferred from summed probability distributions of radiocarbon dates. The Journal of Island and Coastal Archaeology:1–19. doi:10.1080/15564894.2019.1605429. [accessed 2019 Nov 21]. Available at https://www.tandfonline.com/doi/full/10.1080/15564894.2019.1605429.Google Scholar
Brown, WA. 2017. The past and future of growth rate estimation in demographic temporal frequency analysis: Biodemographic interpretability and the ascendance of dynamic growth models. Journal of Archaeological Science 80:96–108. doi: 10.1016/j.jas.2017.02.003. [accessed 2020 Mar 5]. Available at http://www.sciencedirect.com/science/article/pii/S0305440317300328.CrossRefGoogle Scholar
Chapman, J. 2020. Forging identities in the prehistory of Old Europe: dividuals, individuals and communities, 7000–3000 BC. Leiden: Sidestone Press.Google Scholar
Chapman, J, Higham, T, Slavchev, V, Gaydarska, B, Honch, N. 2007. Social context of the emergence, development, and abandonment of the Varna cemetery, Bulgaria. European Journal of Archaeology 9(2):157–181. doi: 10.1177/1461957107086121.Google Scholar
Crema, ER, Bevan, A. 2021. Inference from large sets of radiocarbon dates: software and methods. Radiocarbon 63(1):23–39. doi: 10.1017/RDC.2020.95. [accessed 2021 Aug 25]. Available at https://www.cambridge.org/core/journals/radiocarbon/article/inference-from-large-sets-of-radiocarbon-dates-software-and-methods/DE7B3B8298D6086AEBFE306C714E44BC#.CrossRefGoogle Scholar
Crema, ER, Bevan, A, Shennan, S. 2017. Spatio-temporal approaches to archaeological radiocarbon dates. Journal of Archaeological Science 87:1–9. doi: 10.1016/j.jas.2017.09.007. [accessed 2019 Nov 21]. Available at https://linkinghub.elsevier.com/retrieve/pii/S0305440317301310.CrossRefGoogle Scholar
Crema, ER, Habu, J, Kobayashi, K, Madella, M. 2016. Summed probability distribution of 14C dates suggests regional divergences in the population dynamics of the Jomon Period in eastern Japan. Pinhasi R, editor. PLoS One 11(4):e0154809. doi:10.1371/journal.pone.0154809. Available at https://dx.plos.org/10.1371/journal.pone.0154809.CrossRefGoogle Scholar
Crema, ER, Kobayashi, K. 2020. A multi-proxy inference of Jōmon population dynamics using bayesian phase models, residential data, and summed probability distribution of 14C dates. Journal of Archaeological Science 117:105136. doi: 10.1016/j.jas.2020.105136. [accessed 2020 May 3]. Available at http://www.sciencedirect.com/science/article/pii/S0305440320300583.CrossRefGoogle Scholar
Crema, ER, Shoda, S. 2021. A Bayesian approach for fitting and comparing demographic growth models of radiocarbon dates: a case study on the Jomon-Yayoi transition in Kyushu (Japan). PloS One 16(5):e0251695.CrossRefGoogle Scholar
de Souza, JG, Riris, P. 2021. Delayed demographic transition following the adoption of cultivated plants in the eastern La Plata Basin and Atlantic coast, South America. Journal of Archaeological Science 125:105293. doi: 10.1016/j.jas.2020.105293. Available at https://www.sciencedirect.com/science/article/pii/S0305440320302144.CrossRefGoogle Scholar
Downey, SS, Bocaege, E, Kerig, T, Edinborough, K, Shennan, S. 2014. The neolithic demographic transition in Europe: correlation with juvenility index supports interpretation of the summed calibrated radiocarbon date probability distribution (SCDPD) as a valid demographic proxy. PLoS One 9(8):e105730.CrossRefGoogle ScholarPubMed
Downey, SS, Haas, WR, Shennan, SJ. 2016. European Neolithic societies showed early warning signals of population collapse. Proceedings of the National Academy of Sciences USA 113(35):9751–9756. doi: 10.1073/pnas.1602504113.CrossRefGoogle ScholarPubMed
Edinborough, K, Porčić, M, Martindale, A, Brown, TJ, Supernant, K, Ames, KM. 2017. Radiocarbon test for demographic events in written and oral history. Proceedings of the National Academy of Sciences USA 114(47):12436–12441.CrossRefGoogle ScholarPubMed
Frînculeasa, A. 2016. Nordul Munteniei şi cronologia aspectului cultural Stoicani-Aldeni-stratigrafie, elemente de reper şi date radiocarbon din situl de la Mălăieştii de Jos (jud. Prahova). Buletinul Muzeului Judeţean Teleorman 8:59–107.Google Scholar
Goldberg, A, Mychajliw, AM, Hadly, EA. 2016. Post-invasion demography of prehistoric humans in South America. Nature 532(7598):232–235. doi: 10.1038/nature17176. Available at https://www.nature.com/articles/nature17176.CrossRefGoogle ScholarPubMed
Harper, TK. 2019. Demography and climate in Late Eneolithic Ukraine, Moldova, and Romania: Multiproxy evidence and pollen-based regional corroboration. Journal of Archaeological Science: Reports 23:973–982.Google Scholar
Hengl, T. 2006. Finding the right pixel size. Computers & Geosciences 32(9):1283–1298.CrossRefGoogle Scholar
Hervella, M, Rotea, M, Izagirre, N, Constantinescu, M, Alonso, S, Ioana, M, Lazăr, C, Ridiche, F, Soficaru, AD, Netea, MG. 2015. Ancient DNA from south-east Europe reveals different events during early and middle neolithic influencing the European genetic heritage. PloS One 10(6):e0128810.CrossRefGoogle ScholarPubMed
Lazăr, C, Mărgărit, M, Radu, V. 2018. Between dominant ideologies and techno-economical constraints: Spondylus ornaments from the Balkans in the 5th millennium BC. In: Interchange in pre-and protohistory: case studies in Iberia, Romania, Turkey and Israel. BAR Publishing (BAR International Series 2891). p. 5–22.Google Scholar
Marinescu-Bîlcu, S. 2001. Gumelnița culture. A “lost” civilization: Gumelnița CD-ROM is published with the support of the Ministry of Culture and Religious Affairs, Archaeology Service, Archaeological Publications Programme Bucharest, cIMeC. [accessed 2021 Nov 7]. Available at http://old.cimec.ro/Arheologie/gumelnita/gumelnita_engl/cd/default.htm.Google Scholar
Mathieson, I, Alpaslan-Roodenberg, S, Posth, C, Szécsényi-Nagy, A, Rohland, N, Mallick, S, Olalde, I, Broomandkhoshbacht, N, Candilio, F, Cheronet, O. 2018. The genomic history of southeastern Europe. Nature 555(7695):197–203.CrossRefGoogle ScholarPubMed
McLaughlin, TR. 2019. On applications of space–time modelling with open-source 14C age calibration. Journal of Archaeological Method and Theory 26(2):479–501. doi: 10.1007/s10816-018-9381-3. Available at https://doi.org/10.1007/s10816-018-9381-3.CrossRefGoogle Scholar
Palmisano, A, Bevan, A, Shennan, S. 2017. Comparing archaeological proxies for long-term population patterns: An example from central Italy. Journal of Archaeological Science 87:59–72. doi: 10.1016/j.jas.2017.10.001. [accessed 2020 Aug 7]. Available at http://www.sciencedirect.com/science/article/pii/S0305440317301474.CrossRefGoogle Scholar
Palmisano, A, Woodbridge, J, Roberts, CN, Bevan, A, Fyfe, R, Shennan, S, Cheddadi, R, Greenberg, R, Kaniewski, D, Langgut, D, et al. 2019. Holocene landscape dynamics and long-term population trends in the Levant. The Holocene 29(5):708–727. doi: 10.1177/0959683619826642. [accessed 2019 Nov 21]. Available at http://journals.sagepub.com/doi/10.1177/0959683619826642.CrossRefGoogle Scholar
Pebesma, E, Bivand, R, Rowlingson, B, Gomez-Rubio, V, Hijmans, R, Sumner, M, MacQueen, D, Lemon, J, Lindgren, F, O’Brien, J, et al. 2022. SP: classes and methods for spatial data. Available at https://CRAN.R-project.org/package=sp.Google Scholar
Pebesma, E, Graeler, B. 2021. GSTAT: Spatial and spatio-temporal geostatistical modelling, prediction and simulation. [accessed 2021 Jun 25]. Available at https://CRAN.R-project.org/package=gstat.Google Scholar
Petrova, V. 2016. Varna culture: an autonomous phenomenon or a local version of the Kodzhadermen-Gumelnitsa- Karanovo VІ cultural complex. In: Nokolov V, Schier W, editors. Der Schwarzmeerraum vom neolithicum bis in die früheisenzeit (6000–600 V. CHR.). Leidorf: Rahden/Westf. p. 123–128.Google Scholar
Popovici, DN. 2010. Copper Age traditions north of the Danube River. In: Anthony DW, Chi JY, editors. The Lost World of Old Europe. The Danube Valley, 5000–3500 BC. New York: Princeton University Press. p. 91–111.Google Scholar
Porčić, M, Blagojević, T, Pendić, J, Stefanović, S. 2021. The Neolithic Demographic Transition in the Central Balkans: population dynamics reconstruction based on new radiocarbon evidence. Philosophical Transactions of the Royal Society B: Biological Sciences 376(1816):20190712. doi: 10.1098/rstb.2019.0712. [accessed 2021 Apr 5]. Available at https://royalsocietypublishing.org/doi/10.1098/rstb.2019.0712.CrossRefGoogle ScholarPubMed
Porčić, M, Blagojević, T, Stefanović, S. 2016. Demography of the early Neolithic population in central Balkans: population dynamics reconstruction using summed radiocarbon probability distributions. PLoS One 11(8):e0160832.CrossRefGoogle ScholarPubMed
Porčićić, M, Nikolić, M. 2016. The approximate Bayesian computation approach to reconstructing population dynamics and size from settlement data: demography of the Mesolithic-Neolithic transition at Lepenski Vir. Archaeological and Anthropological Sciences 8(1):169–186. doi: 10.1007/s12520–014-0223-2. [accessed 2021 Jul 23]. Available at https://doi.org/10.1007/s12520–014-0223-2.CrossRefGoogle Scholar
R Core Team. 2020. R: A language and environment for statistical computing. R Foundation for Statistical Computing. Available at http://www.R-project.org/.Google Scholar
Reimer, PJ, Austin, WE, Bard, E, Bayliss, A, Blackwell, PG, Ramsey, CB, Butzin, M, Cheng, H, Edwards, RL, Friedrich, M, et al. 2020. The IntCal20 Northern Hemisphere radiocarbon age calibration curve (0–55 cal kBP). Radiocarbon 62(4)725–757.CrossRefGoogle Scholar
Reingruber, A, Thissen, L. 2017. The 14SEA Project. A 14C database for southeast Europe and Anatolia (10,000–3000 calBC). [accessed 2019 Sep 8]. Available at http://www.14sea.org/index.html.Google Scholar
Riris, P. 2018. Dates as data revisited: A statistical examination of the Peruvian preceramic radiocarbon record. Journal of Archaeological Science 97:67–76. doi: 10.1016/j.jas.2018.06.008. [accessed 2019 Nov 21]. Available at https://linkinghub.elsevier.com/retrieve/pii/S0305440318301067.CrossRefGoogle Scholar
Roberts, N, Woodbridge, J, Bevan, A, Palmisano, A, Shennan, S, Asouti, E. 2018. Human responses and non-responses to climatic variations during the last Glacial-Interglacial transition in the eastern Mediterranean. Quaternary Science Reviews 184:47–67.CrossRefGoogle Scholar
Sakamoto, Y, Ishiguro, M, Kitagawa, G. 1986. Akaike information criterion statistics. Dordrecht, Netherlands: D. Reidel.Google Scholar
Shennan, S, Downey, SS, Timpson, A, Edinborough, K, Colledge, S, Kerig, T, Manning, K, Thomas, MG. 2013. Regional population collapse followed initial agriculture booms in mid-Holocene Europe. Nature Communications 4. doi:10.1038/ncomms3486. [accessed 2014 Aug 23]. Available at http://www.nature.com/doifinder/10.1038/ncomms3486.CrossRefGoogle Scholar
Slavchev, V. 2010. The Varna Eneolithic cemetery in the context of the late copper age in the east Balkans. In: Anthony DW, Chi J, editors. The lost world of Old Europe. The Danube Valley, 5000–3500 BC. New York: Institute for the Study of the Ancient World. p. 192–211.Google Scholar
Sterud, E, Evans, RK, Rasson, JA. 1984. Part I: Ex Balcanis Lux? Recent developments in Neolithic and Chalcolithic research in southeast Europe. American Antiquity. p. 713–741.Google Scholar
Timpson, A, Barberena, R, Thomas, MG, Méndez, C, Manning, K. 2020. Directly modelling population dynamics in the South American arid diagonal using 14C dates. Philosophical Transactions of the Royal Society B: Biological Sciences 376(1816):20190723. doi: 10.1098/rstb.2019.0723.CrossRefGoogle Scholar
Timpson, A, Colledge, S, Crema, E, Edinborough, K, Kerig, T, Manning, K, Thomas, MG, Shennan, S. 2014. Reconstructing regional population fluctuations in the European Neolithic using radiocarbon dates: a new case-study using an improved method. Journal of Archaeological Science 52:549–557.CrossRefGoogle Scholar
Timpson, A, Manning, K, Shennan, S. 2015. Inferential mistakes in population proxies: A response to Torfing’s “Neolithic population and summed probability distribution of 14C-dates.” Journal of Archaeological Science 63:199–202. doi: 10.1016/j.jas.2015.08.018.CrossRefGoogle Scholar
Todorova, H. 1986. Kamenno-mednata Epokha v Bulgariya. Peto Khilyadoletie Predi Novata Era. Sofia: Nauka i Izkustvo.Google Scholar
Todorova, H. 2003. Praehistory of Bulgaria. In: Grammenos, DV, editor. Recent research in the prehistory of the Balkans. Thessaloniki: Archaeological Institute of Northern Greece. p. 257–328.Google Scholar
Todorova, K, Zhelyaskova, V. 1978. The Eneolithic period in Bulgaria in the fifth millennium B.C. Oxford: BAR (BAR International Series). [accessed 2021 Nov 13]. Available at https://ehrafarchaeology.yale.edu/document?id=e075-006.Google Scholar
Vander Linden, M, Silva, F. 2021. Dispersals as demographic processes: testing and describing the spread of the Neolithic in the Balkans. Philosophical Transactions of the Royal Society B 376(1816):20200231.CrossRefGoogle ScholarPubMed
Vrhovnik, DM. 2019. Neolithic and Copper Age settlement dynamics in the Western Carpathian Basin and Eastern Alps. Documenta Praehistorica 46:268–282.CrossRefGoogle Scholar
Weninger, B, Clare, L, Jöris, O, Jung, R, Edinborough, K. 2015. Quantum theory of radiocarbon calibration. World Archaeology 47(4):543–566. doi: 10.1080/00438243.2015.1064022. [accessed 2020 Sep 23]. Available at https://doi.org/10.1080/00438243.2015.1064022.CrossRefGoogle Scholar
Weninger, B, Joris, O, Danzeglocke, U. 2020. CalPal Holocene palaeolithic 14C database. www.calpal.de. [accessed 2021 Apr 13]. Available at https://www.academia.edu/40774947/CalPal_Holocene_Palaeolithic_14C_Database.Google Scholar
Figure 0

Figure 1 Map of radiocarbon distribution data set. Legend: 1-Akladi Cheiri; 2-Bikovo; 3-Čardako-Slatino; 4-Djakovo; 5-Dolnoslav; 6-Drama-Merdžumekja; 7-Durankulak; 8-Ezero; 9-Goljamo Delčevo; 10-Hotnica; 11-Junacite; 12-Karnobat; 13-Košarna; 14-Omurtag; 15-Orlitsa; 16-Ovčarovo; 17-Povelyanovo; 18-Smjadovo; 19-Sušina; 20-Tatul; 21-Tell Azmak; 22-Tell Karanovo; 23-Tell Russe; 24-Varhari; 25-Varna1; 26-Varna2; 27-Varna3; 28-Baia Boruz Tell; 29-Popina Blagodeasca; 30-Bordușani; 31-Carcaliu; 32-CăscioareleOstrovel; 33-Cunești; 34-Dambul lui Haralambie; 35-Gumelnița-terrasse; 36-Gumelnița-tell; 37-Hârșova; 38-Lișcoteanca-Movila Olarului; 39-Lunca; 40-Luncavița; 41-Mălăieștii de Jos; 42-Măriuța-C; 43-Măriuța-T; 44-Navodari; 45-Niculițel; 46-Orbeasca Sus; 47-Panduru; 48-Pietrele; 49-Seciu; 50-Șeinoiu; 51-SultanaGhețărie; 52-Sultana-Malu Roșu-terrasse; 53-Sultana-Malu Roșu-tell; 54-Taraschina; 55-Taraschina_2; 56-Urlați; 57-Vărăști; 58-Vitănești; 59-Vlădiceasca.

Figure 1

Figure 2 (a) Number of radiocarbon dates per grid cell. Values are log10 scaled. (b) The earliest appearance of the KGK-VI settlements in the grid cells, consisting of grid cell centroids with the date for the beginning of the KGK-VI occupation. Gridded area (white hexagons) represents the KGK-VI area with dated sites.

Figure 2

Figure 3 Results of fitting and comparing the entire regional empirical SPD against the exponential null model of population growth. Monte Carlo 95% confidence null model gray envelope is based on 5000 runs. Observed SPD is shown with solid red line, while the positive and negative deviations from the null are marked in red and blue. (Please see online version for color figures.)

Figure 3

Figure 4 Results of fitting and comparing the entire regional empirical SPD against the logistic null model of population growth. Monte Carlo 95% confidence null model gray envelope is based on 5000 runs. Observed SPD is shown with solid red line, while the positive and negative deviations from the null are marked in red and blue.

Figure 4

Figure 5 Bootstrapped composite kernel density estimate, suggesting a composite model with the breakpoint at approximately 4400 BC should be tested.

Figure 5

Figure 6 The KGK-VI empirical SPD record fitted to logistic (5000 BC–4400 BC) and exponential (4400 BC–5750 BC) models, with significance envelope derived from 5000 Monte Carlo simulations.

Figure 6

Figure 7 Permutation test showing variation between regional population growth. Observed SPDs for each region are shown with a solid black line, while the dashed line represents the observed pan-regional SPD. Gray areas represent the 95% confidence envelope for the null model, red and blue bands represent areas where the observed SPD significantly positively (red) and negatively (blue) deviates from the pan regional null model.

Figure 7

Figure 8 Observed rate of growth at each transition computed from the SPD. I: 5050–4800 to 4800–4550; II: 4800–4550 to 4550–4300; III: 4550–4300 to 4300–4050; IV: 4300–4050 to 4050–3800 BC.

Figure 8

Figure 9 Local geometric growth rate for each transition block. I: 5050–4800 to 4800–4550; II: 4800–4550 to 4550–4300; III: 4550–4300 to 4300–4050; IV: 4300–4050 to 4050–3800 BC; V: Geographical reference map for the local geometric growth rate analysis, shown in transitional blocks I-IV.

Figure 9

Figure 10 Spatial permutation test showing where growth is significantly higher or lower than the null for each transition block. I:5050–4800 to 4800–4550; II:4800–4550 to 4550–4300; III:4550–4300 to 4300–4050; IV:4300–4050 to 4050–3800 BC; V: Geographical reference map for the spatial permutation test analysis, shown in transitional blocks I-IV. Significance is shown in terms of q-values (more robust against false positives) and p-values.

Supplementary material: File

Popescu et al. supplementary material

Popescu et al. supplementary material

Download Popescu et al. supplementary material(File)
File 483.8 MB