Ecosystem productivity affected the spatiotemporal disappearance of Neanderthals in Iberia
Fauna, culture and chronology datasetsA geo-referenced dataset of chronometric dates covering the late MIS 3 (55–30 kyr cal bp) was compiled from the literature (dataset 1). The dataset included 363 radiocarbon, thermoluminescence, optically stimulated luminescence and uranium series dates obtained from 62 archaeological sites and seven palaeontological sites. These chronological determinations were obtained from ten palaeontological levels and 138 archaeological levels. The archaeological levels were culturally attributed to the Mousterian (n = 75), Châtelperronian (n = 6) and Aurignacian (n = 57) technocomplexes. A number of issues can potentially hamper the chronological assessment of Palaeolithic technocomplexes from radiocarbon dates, such as pretreatment protocols that do not remove sufficient contaminants or the quality of the bone collagen extracted. Moreover, discrepancies in cultural attributions or stratigraphic inconsistencies are commonly detected in Palaeolithic archaeology. Information regarding the quality of date determinations and cultural attribution or stratigraphic issues is provided in the Supplementary Information.Our dataset also included the presence of herbivore species recovered from each archaeo-palaeontological site (hereafter referred to as local faunal assemblages (LFAs)), their body masses and their chronology. The mean body mass of both sexes, for each species, was obtained from the PHYLACINE database53 and used in the macroecological modelling approach described below (see ‘Carrying capacity of herbivores’). For visual representation purposes, the herbivore species were grouped into four weight categories: small (500 kg). The chronology of the occurrence of each herbivore species was assumed to be the same as the dated archaeo-palaeontological layer where the species remains were recovered. Thus, to estimate the chronological range of each species in each region, all radiocarbon determinations were calibrated with the IntCal20 calibration curve54 and OxCAL4.2 software55. The BAMs were run to compute the upper and lower chronological boundaries at a CI of 95.4% of each LFA (see ‘Chronological assessment’ for more details). One of the purposes of the current study was to estimate the potential fluctuations in herbivore biomass during the stadial and interstadial periods of the late MIS 3. Accordingly, the time spans of the LFAs were classified into the discrete GS and GI phases provided by Rasmussen et al.51.Geographic settingsThe Iberian Peninsula locates at the southwestern edge of Europe (Fig. 1). It constitutes a large geographic area that exhibits a remarkable diversity of ecosystems, climates and landscapes. Both now and in the past, altitudinal, latitudinal and oceanic gradients affected the conformation of two biogeographical macroregions with different flora and fauna species pools: the Eurosiberian and Mediterranean regions13,46. In the north, along the Pyrenees and Cantabrian strip, the Eurosiberian region is characterized by oceanic influence and mild temperatures in the present day, whereas the Mediterranean region features drier summers and milder winters (Fig. 1). Between the Eurosiberian and Mediterranean regions, there is a transitional area termed Submediterranean or Supramediterranean. Lastly, the Mediterranean region is divided into two distinctive bioclimatic belts: (1) the Thermomediterranean region, located at lower latitudes, with high evapotranspiration rates and affected by its proximity to the coast; and (2) the Mesomediterranean region, with lower temperatures and wetter conditions (Fig. 1).Previous studies have shown that zoocoenosis and phytocenosis differed between these macroregions in the Pleistocene13,46. However, flora and fauna distributions changed during the stadial–interstadial cycles in the Iberian Peninsula, which suggests potential alterations in the boundaries of these biogeographical regions. The modelling approach used in this study to estimate the biomass of primary consumers is dependent on the reconstructed NPP and the herbivore guild structure in each biogeographical region. To test the suitability of the present-day biogeographical demarcations of the Iberian Peninsula during MIS 3, we assessed whether the temporal trends of NPP and the composition of each herbivore palaeocommunity differed between these biogeographical regions during the MUPT.Chouakria and Nagabhusan56 proposed a dissimilarity index to compare time series data by taking into consideration the proximity of values and the temporal correlation of the time series:$${rm{CORT}}(S_1,S_2) = frac{{mathop {sum}nolimits_{i = 1}^{p – 1} {left( {u_{left( {i + 1} right)} – u_i} right)} (v_{(i + 1)} – v_i)}}{{sqrt {mathop {sum}nolimits_{i = 1}^{p – 1} {(u_{(i + 1)} – u_i)^2} } sqrt {mathop {sum}nolimits_{i = 1}^{p – 1} {(v_{(i + 1)} – v)^2} } }}$$
(1)
where S1 and S2 are the time series of data, u and v represent the values of S1 and S2, respectively, and p is the length of values of each time series. CORT(S1, S2) belongs to the interval (−1,1). The value CORT(S1, S2) = 1 indicates that in any observed period (ti, ti+1), the values of the sequence S1 and those of S2 increase or decrease at the same rate, whereas CORT = −1 indicates that when S1 increases, S2 decreases or vice versa. Lastly, CORT(S1, S2) = 0 indicates that the observed trends in S1 are independent of those observed in S2. To complement this approach by considering not only the temporal correlation between each pair of time series but also the proximity between the raw values, these authors proposed an adaptive tuning function defined as follows:$$d{rm{CORT}}left( {S_1,S_2} right) = fleft({{rm{CORT}}left( {S_1,S_2} right)} right)times dleft( {S_1,S_2} right)$$
(2)
where$$fleft( x right) = frac{2}{{1 + exp left( {k,x} right)}},k ge 0$$
(3)
In this study, k was 2, meaning that the behaviour contribution was 76% and the contribution of the proximity between values was 24%57. Hence, f(x) modulates a conventional pairwise raw data distance (d(S1,S2)) according to the observed temporal correlation56. Consequently, dCORT adjusts the degree of similarity between each pair of observations according to the temporal correlation and the proximity between values. This function was used to compare the reconstructed NPP between biogeographical regions during MIS 3 in the Iberian Peninsula. However, two different biogeographical regions could have experienced similar evolutionary trends in their NPP, even though their biota composition was different. Therefore, this analysis was complemented with a JSI to assess whether the reconstructed herbivore species composition in each palaeocommunity differed among biogeographical regions during the late MIS 3. The JSI was based on presence–absence data and was calculated as follows:$${rm{JSI}} = frac{c}{{(a + b + c)}}$$
(4)
where c is the number of shared species in both regions and a and b are the numbers of species that were only present in one of the biogeographical regions. Therefore, the higher the value the more similar the palaeocommunities of both regions were.Chronological assessmentPivotal to any hypothesis of Neanderthal replacement patterns by AMHs is the chronology of that population turnover. To this end, we used three different approaches to provide greater confidence in the results: BAMs, the OLE model and SPD of archaeological assemblages. As detailed below, each of these approaches provides complementary information about the MUPT.First, we built a set of BAMs for the Mousterian, Châtelperronian and Aurignacian technocomplexes in each region during the MIS 3. As stated above, we compiled the available radiocarbon dates for Iberia between 55 and 30 kyr cal bp. However, not all dates or levels were included in the Bayesian chronology models. Radiocarbon determinations obtained from shell remains were incorporated in the dataset (dataset 1); however, the local variation of the reservoir age was unknown from 55 to 30 kyr bp. Because of uncertainties related to marine reservoir offsets, all BAMs that incorporated dates from marine shells were run twice: including and excluding these dates. All of the archaeological levels with cultural attribution issues or stratigraphic inconsistencies were excluded. The Supplementary Note provides a detailed description of the sites, levels and dates excluded and their justification. All BAMs were built for each technocomplex using the OxCAL4.2 software55 and IntCal20 calibration curve54.Bayesian chronology models were built for each archaeological and palaeontological level. Then, the dates associated with each technocomplex were grouped within a single phase to determine each culture’s regional appearance or disappearance. Our interest was not focused on the chronological duration of the Mousterian, Châtelperronian and Aurignacian cultures, but on the probability distribution function of the temporal boundaries of these cultures in each region. Thus, this chronological assessment aims to provide an updated chronological frame for Neanderthal replacement by AMHs in Iberia. For this reason, we did not differentiate between proto- and early Aurignacian cultures, since both are attributed to AMHs.In each BAM, we inserted into the same sequence the radiocarbon dates associated with a given technocomplex within a start and end boundary to bracket each culture, which allowed us to determine the probability distribution function for the beginning and end moment of each cultural phase6. The resolution of all models was set at 20 years. We used a t-type outlier model with an initial 5% probability for each determination, but when more than one radiocarbon date was obtained from the same bone remain, we used an s-type outlier model and the combine function. The thermoluminescence dating likelihoods were included in the models, together with their associated 1σ uncertainty ranges. When dates with low agreement ( More
