Abstract The replacement of Neanderthals by Anatomically Modern Humans has typically been attributed to environmental pressure or a superiority of modern humans with respect to competition for resources. Here we present two independent models that suggest that no such heatedly debated factors might be needed to account for the demise of Neanderthals. Starting from the observation that Neanderthal populations already were small before the arrival of modern humans, the models implement three factors that conservation biology identifies as critical for a small population’s persistence, namely inbreeding, Allee effects and stochasticity. Our results indicate that the disappearance of Neanderthals might have resided in the smallness of their population(s) alone: even if they had been identical to modern humans in their cognitive, social and cultural traits, and even in the absence of inter-specific competition, Neanderthals faced a considerable risk of extinction. Furthermore, we suggest that if modern humans contributed to the demise of Neanderthals, that contribution might have had nothing to do with resource competition, but rather with how the incoming populations geographically restructured the resident populations, in a way that reinforced Allee effects, and the effects of inbreeding and stochasticity.

Citation: Vaesen K, Scherjon F, Hemerik L, Verpoorte A (2019) Inbreeding, Allee effects and stochasticity might be sufficient to account for Neanderthal extinction. PLoS ONE 14(11): e0225117. https://doi.org/10.1371/journal.pone.0225117 Editor: Sergi Lozano, Universitat de Barcelona, SPAIN Received: March 18, 2019; Accepted: October 29, 2019; Published: November 27, 2019 Copyright: © 2019 Vaesen et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: Data will be made available at DOI: 10.17026/dans-znu-2ssx. Funding: Research by KV was supported by the Netherlands Organisation for Scientific Research (NWO) (grant 276-20-021). Competing interests: The authors have declared that no competing interests exist.

Introduction A long-standing enigma in palaeoanthropology is the demise of Neanderthals approximately 40 kya [1]. There is general agreement that―after a long period of largely separated coexistence [2–6]―their disappearance roughly coincides with migration events starting around ~60 kya by Anatomically Modern Humans (AMHs) from Africa into the Near East and Europe and that, accordingly, AMHs took over the territories previously occupied by our sister species [7–12]. What is uncertain however, are the causes of Neanderthal extinction. It has been attributed to a wide variety of intensely debated factors, including climatic change ([13–17], but see [18–19]), epidemics [20–21], a superiority of AMHs over Neanderthals in competing for the same resources ([22–34] but see [18]). Models have, not surprisingly, confirmed that if such superiority is assumed [35–41], Neanderthals would indeed have been likely to go extinct. Here we argue that no such contested factors might be needed to account for the demise of Neanderthals. We present two independent models that capture the internal dynamics of Neanderthal populations―the models thus ignore, among other things, competitive interactions with AMHs―and that suggest that the disappearance of Neanderthals might have resided in the small size of their population(s) alone. Accordingly, our study substantiates the suggestion, made in passing by French [42], that “it may simply be the case that Neanderthal populations declined below their minimum viable population threshold”. Genetic studies, indeed, indicate that Neanderthal effective population size―the size of the ideal population that would undergo the same amount of random genetic drift as the actual population [43]―was already small by ~400 kya, amounting to no more than 3,000–3,900 individuals, a level that was sustained almost up till the species’ extinction ~38 kya [44–46]. Rogers et al. [47] report a higher number for the effective size of the global Neanderthal population, but also show that this global population was subdivided in very small and highly isolated local populations (each, so to speak, with its own, small effective size). Likewise, Castellano et al. [48] conclude, based on a comparison between the genome of three Neanderthal individuals and the genome of present-day modern humans, that Neanderthal genetic diversity was limited and that Neanderthal populations were small and isolated from one another. The archaeological data point in the same direction ([49], and references therein); and Bocquet-Appel and Degioanni [50] estimate that the Neanderthal census metapopulation was in the range of a mere 5,000–70,000 individuals. The models presented here implement three basic factors that, according to conservation biology (the field from which our models were drawn), would put such small populations at risk of extinction: inbreeding, Allee effects and stochasticity ([51]; in archaeology, see Finlayson [52]). Inbreeding depression refers to the reduction in fitness of individuals that arise from matings between genetic relatives, matings thus that are more likely to occur in small populations. Inbreeding, which seems to have been common in Neanderthals [44, 53–55], might lead to a lower fitness because it increases the chances of the expression of recessive, deleterious traits and because homozygotes often have a general disadvantage relative to heterozygotes. Harris and Nielsen [56] estimate that, due to inbreeding, Neanderthals had at least 40% lower fitness than modern humans on average. Allee effects refer to the effects that population density has on reproduction and, thus, on population growth [51]. At lower densities, the case we are concerned with here, growth rates might drop due to problems in mate-finding, and to several problems that highly cooperative species, such as Neanderthals, are particularly susceptible to, including low availability of helpers in cooperative hunting, defending kills from kleptoparasites, and allo-parenting [57]. Finally, stochastic, annual fluctuations in births, deaths and sex ratio are more likely to place smaller populations on a trajectory towards extinction than bigger ones [51]. Our models indicate that these factors alone could have resulted in Neanderthal extinction, even if Neanderthals and AMHs were identical in terms of individual-level traits that are deemed relevant to persistence or extinction (e.g., cognitive and technological ability, sociality).

Materials and methods The specific question our study aims to address is whether inbreeding, Allee effects and stochasticity are sufficient to explain the disappearance of Neanderthals. To that end, we develop two separate models, both of which track Neanderthal population growth over time. Crucially, since we want to avoid making assumptions about the superiority of AMHs over Neanderthals, both models are parameterized based on estimates for AMHs, estimates that pertain to, among other things, reproduction, mortality, and inbreeding (see below, S1 Appendix, and S6 and S7 Tables). The first model is a deterministic matrix model, the second an individual-based stochastic simulation model (IBM). The matrix model served two purposes: to calibrate the IBM and to validate some of the results of the IBM. Since, in contrast to the matrix model, the IBM allowed us straightforwardly to introduce inbreeding and stochasticity, it was our primary resource for obtaining the results reported here. More specifically, we used the IBM to determine, separately, the levels of inbreeding and the levels of Allee effects that put Neanderthal populations at risk of extinction; and to determine Neanderthal’s vulnerability to extinction in a scenario involving both inbreeding and Allee effects. Subsequently, we set these results against what―as it comes to Neanderthal population size, inbreeding and Allee effects―can be inferred from the literature. Below we provide a general, non-technical overview of the research set-up. For a full description, please consult S1 Appendix. General description of the basic models The deterministic model consists of a Leslie matrix, a type of matrix that is commonly used in conservation biology to model population growth [58]. The matrix summarizes, for each of a population’s age classes, yearly survival and reproduction. Our model assumes a sex ratio of 1:1 and, because they are the limiting factor for reproduction, only considers females. Accordingly, the first row of the matrix contains, for each of the age classes, the average number of female offspring per female per year. The subdiagonals represent the yearly survivals of females; these values are on the subdiagonals because, whenever a female survives in a given year, she will move to the next age class. The growth factor of the population―the factor by which the population will annually be multiplied―is given by the dominant eigenvalue of the matrix. If this factor is less than 1, the population will have a negative growth rate and goes extinct. For the IBM, we relied on VORTEX [59–61], a software package used by conservation biologists to perform Population Viability Analyses of endangered wildlife species. VORTEX simulates the annual life events (e.g., sex determination, breeding, mortality) that might occur to each of the individuals within a given population, and records over time, in discrete time steps, the characteristics of these individuals as well as those of the population as a whole. Occurrences of events are probabilistic. Demographic stochasticity (relating to annual fluctuations in, e.g., sex ratios, births and deaths) is thus inherent to VORTEX models. Calibration We calibrated the two models by checking whether, under similar parameter settings, they were able to produce similar trends in Neanderthal population growth. At this stage, we didn’t include inbreeding or Allee effects. The primary input parameters for both models were derived from estimates of female reproduction in extant hunter-gatherers [62] and from the West model life table, Level 5 [63]. In general, such model life tables summarize the mortality events of an ideal population. West tables are the most commonly used in human palaeo-demography; Level 5 is closely matched by the mortality profiles of extant hunter-gatherers [64]. The parameters used in the matrix model are shown in S1 Appendix, S1 Table. Calculating the dominant eigenvalue of the corresponding matrix yields a population growth factor of 1.008, which corresponds to a relative growth rate of 0.80%. The parameters used in the VORTEX model are summarized in S1 Appendix, S5 and S6 Tables. For various initial population sizes (N 0 = 50; 100; 500; 1,000; 5,000), we ran ten simulation runs, each simulating a time span of 100 years. The relative growth rate of these simulations was, on average, 0.76%. This value is slightly lower than the relative growth rate obtained from the matrix model. However, it has been theoretically established that, in general, growth rates of stochastic models are less than the growth rates of deterministic models [65]. Inbreeding and stochasticity In conservation biology there is no standard procedure for implementing inbreeding in deterministic matrix models. Therefore, we relied on our IBM to model inbreeding. In VORTEX, inbreeding depression is modeled in terms of its effects on infant survival. It is governed by two parameters: the number of lethal equivalents, I, which is the number of recessive alleles carried in a heterozygous genome that would be lethal if carried in the homozygous state; and the percentage, f i , of the inbreeding depression caused by such lethal alleles rather than by other genetic mechanisms (e.g., a general disadvantage of homozygotes). Based on observations of several species, VORTEX sets f i at 50% by default. For various initial population sizes (N 0 = 50; 100; 500; 1,000; 5,000) and various value of f i (viz., f i = 30;50;70), and again starting from West table Level 5, we determined, by simulating widely over I, the parameter I risk . I risk is the lowest value of I that yields at least one extinction event in ten simulation runs, each run simulating over a time span of 10,000 years. Furthermore, we determined, I sure , which is the lowest value of I that yields extinction in all of the ten runs. Note that VORTEX models that incorporate inbreeding generally run very slow, and do so especially when carrying capacity is high. Therefore we set K at 5,000, unlike we did in the basic model (where carrying capacity was set at K = 10,000). This doesn’t affect our conclusions (see below). Allee effects, with and without stochasticity The dependence of female breeding rates on population size is commonly (including in VORTEX) captured by an equation that, when simplified to cover the most conservative scenario, tells us that the fraction of females breeding at population size N, E, is given by (1) where p 0 is the percentage of females breeding in the absence of population size constraints, and A is a parameter describing the strength of the Allee effect (see equations S4 and S5 in the S1 Appendix). Accordingly, in the IBM we simulated widely over A, and determined, for various initial population sizes (N 0 = 50; 100; 500; 1,000; 5,000), A risk , which is the lowest value of A that leads, in ten simulation runs (each run comprises 10,000 years), to at least one extinction event. Also, we determined A sure , which is the lowest value of A for which all ten runs result in extinction. In order to validate our results, we introduced the Allee equation used in the IBM into the matrix model, and set the results against the results obtained from the IBM. Inbreeding, Allee effects and stochasticity Estimates of I for modern humans [66–70] range from 0.58 [69] to 2.2 [66]. Gao et al. [69] point out that, since these estimates are based on reported deaths after birth and thus do not take into account prenatal deaths, the actual number of lethal equivalents might be higher; the authors surmise that prenatal deaths might increase I by one additional lethal equivalent (resulting in a maximum value of I = 3.2). To assess the combined effects of inbreeding, Allee effects and stochasticity, we ran simulations for the highest I-value just reported (i.e., I = 3.2). More specifically, by widely varying over A, we determined for this I-value, and for various initial population sizes (N 0 = 50; 100; 500), A risk . Although I = 3.2 is the least conservative value among the values found in the literature, our choice was motivated by the results obtained in scenarios that involved inbreeding alone (i.e., the results suggested that even at I = 3.2 the impact of inbreeding would be relatively small; see Results), and by the large computational demands of VORTEX in scenarios that combine Allee effects, inbreeding and very lengthy timespans (10,000 years). Also due to resource constraints, we did not determine A sure , restricted N to the range 0–500, and set carrying capacity at K = 5,000.

Conclusion and discussion Our results support the hypothesis that the disappearance of Neanderthals might have been the result of a demographic factors alone, that is, the result merely of the internal dynamics that operate in small populations. Our conclusions are consistent with but go beyond the conclusions of a recent study by Kolodny and Feldman [73]. Based on a series of mathematical models, these authors too argue that no external factors (climate, epidemics) nor a superiority of AMHs in resource competition are needed to account for Neanderthal extinction. Their models suggest that migratory dynamics―with more migration happening from Africa into Europe by AMHs than migration from Europe into Africa by Neanderthals―might have been sufficient to result in the replacement of Neanderthals by AMHs. While Kolodny and Feldman’s models indeed do not assume a competitive advantage for either species (but see S1 Appendix), they do take for granted that Neanderthals and AMHs competed for the same habitats. Our study shows that even without this contested assumption ([74]; see also the literature on competition avoidance among extant hunter-gatherers, e.g., [75] and references therein), Neanderthal extinction might have taken place. If Neanderthals lived in small populations since ~400 kya [44–46], why did it take so long for them to become extinct? A first relevant consideration concerns demographic stochasticity. We have seen that annual fluctuations in births, deaths and sex ratio might determine whether and when a small population disappears. So our results are consistent with a scenario in which a small population of Neanderthals persists for several thousands of years, and then, due to a stroke of bad luck, disappears. Furthermore, for the sake of simplicity, our models do not take into account environmental stochasticity. That is, the models work with fixed probabilities for mortality, fertility and sex ratios―fluctuations are thus simply caused by probabilistic sampling (demographic stochasticity). In natural conditions, though, the probabilities themselves will vary, according to random fluctuations in the environment (environmental stochasticity). Note that these fluctuations do not correspond to the millennial trends observed in the palaeoclimatic record, but occur at much lower temporal scales (e.g., a couple of years of drought, an epidemic among prey). They should be understood as natural variations around a mean, rather than as external forcings, such as the ones that some scholars have claimed to be responsible for the demise of Neanderthals (e.g., climatic change [13–16] or volcanic eruptions [17]). Importantly, in a given year, demographic and environmental stochasticity might very well work in opposite directions; environmental conditions, for instance, might be favorable and alleviate the stress induced by demographic stochasticity. In fact, it will very rarely happen that a metapopulation comprising several sub-populations has all the stochastic odds against it, that is, it will only very rarely happen that, for a significant amount of time, environmental variability produces low fertility rates and high death rates, and additionally, demographic stochasticity produces low fertility rates and high death rates and unfavorable sex ratios, and this in all of the metapopulation’s sub-populations. But in the very long run, such an unfavorable scenario eventually will take place. Accordingly, it is not implausible that, despite regular local extinction events [76], a small metapopulation manages to survive over prolonged stretches of time but eventually dies out due to its overall size and stochasticity. Noteworthy, there is nothing unusual about the persistently small size of Neanderthal populations. Hominin populations likely were small throughout the Pleistocene [77]. We suggest that AMHs might still have contributed to the extinction, but not necessarily by engaging in competition with or outcompeting Neanderthals. The mere interspersal of AMH sub-populations between Neanderthal sub-populations reduced the opportunities for intra-breeding and migratory activity among the latter. The resulting small and isolated sub-populations (as documented by [47–49,76,78]) would be increasingly vulnerable precisely to the factors examined in the current paper (viz. inbreeding, Allee effects and stochasticity), and thus to extinction. As such, the presence of modern humans in Eurasia would have accelerated a process that, at some point, was likely to have occurred anyway. Stated otherwise, the arrival of AMHs would have been a contributory factor rather than the cause of the extinction. Importantly, population-level characteristics―e.g., many of the characteristics that conservation biology has shown to be critical for a species’ persistence, including population size, species distribution, intraspecific variability, and patterns of dispersal―might also account for the successful range expansion of AMHs. In other words, our species’ success need not be the result of a superiority in its individual-level traits. An explanation solely in terms of the internal dynamics of the Neanderthal population, as the one presented here, serves as a null hypothesis against which competing, and less parsimonious, hypotheses are to be assessed. Regardless of whether external factors (climate or epidemics) or factors related to resource competition played a role in the actual demise of Neanderthals, our study suggests that any plausible explanation of the demise also needs to incorporate demographic factors as key variables.

Acknowledgments We thank the members of the Human Origins Group, Faculty of Archaeology, Leiden University, and two anonymous reviewers for their feedback on earlier drafts of this paper.