Jhelam N Deshpande

Oliver Kaltz

Emanuel A Fronhofer
email: emanuel.fronhofer@umontpellier.fr

Host-parasite dynamics set the ecological theatre for the evolution of state-and context-dependent dispersal in hosts

Keywords: dispersal evolution, information use, infection, epidemiology, condition-dependent dispersal, migration

While host-parasite interactions are ubiquitous, the large scale consequences of parasite infections are mainly driven by the spatial context. One trait of pivotal importance for the eco-evolutionary dynamics of such metapopulations is the spatial behaviour of hosts, that is, their dispersal. It is well established that dispersal is not a random process, rather dispersal is informed and may depend on internal and external factors. In host-parasite metapopulations, dispersal may be a function of a host's infection state, but also of the local context, such as host density or parasite prevalence. Using a dynamical host-parasite metapopulation model of a parasite that reduces host fecundity (sterilizing parasite), we explore whether host dispersal evolves to be state-and context-dependent and what shapes the evolutionarily stable dispersal reaction norms have. We show that state-dependent dispersal readily evolves in the sense that hosts disperse more when infected, except when infected hosts pay significantly higher dispersal costs. This dispersal bias evolves due to kin selection which is consistent with previous studies. Most importantly, we show that prevalence-dependent dispersal evolves, especially when virulence is high and epidemiological dynamics have predictable signatures. The observed evolutionary outcome, a negatively prevalencedependent dispersal reaction norm for susceptible hosts at high virulence, seems counter-intuitive at first. However, our results can be readily explained by the emergent epidemiological dynamics, especially their spatial and temporal correlation patterns. Finally, we show that context-dependency in dispersal may rely on both, prevalence, but also host density cues. Our work provides new insights into the evolution of complex dispersal phenotypes in host-parasite metapopulations as well as on associated consequences of ecological dynamics for evolutionary change.

Introduction

In many host-parasite systems, spatial behaviour of hosts, more precisely, their movement and dispersal, is key to the spread of infection, especially if populations are spatially structured. In such spatially structured populations (metapopulation sensu lato; see, e.g., [START_REF] Hanski | Metapopulation Ecology[END_REF], the landscape is simplified and represented as a network of interconnected patches (nodes), where the probability of local transmission depends, for instance, on the global rates at which infected and uninfected individuals move between patches, the direction they choose, or the distances they cover (e.g., [START_REF] Kerr | Local migration promotes competitive restraint in a host-pathogen 'tragedy of the commons[END_REF][START_REF] Ferrari | The dynamics of measles in sub-Saharan Africa[END_REF][START_REF] Metcalf | Persistence in epidemic metapopulations: quantifying the rescue effects for measles, mumps, rubella and whooping cough[END_REF][START_REF] Yashima | Epidemic process over the commute network in a metropolitan area[END_REF][START_REF] Dudas | Virus genomes reveal factors that spread and sustained the Ebola epidemic[END_REF]. A metapopulation approach of disease transmission is not only important for understanding and predicting the spread of epidemics across landscapes [START_REF] Ostfeld | Spatial epidemiology: an emerging (or re-emerging) discipline[END_REF][START_REF] Keeling | Networks and epidemic models[END_REF][START_REF] North | The dynamics of disease in a metapopulation: The role of dispersal range[END_REF]. Variation in host and parasite dispersal may also have a strong impact on the evolution of attack and defence traits, such as virulence and resistance, or patterns of local adaptation (reviewed in [START_REF] Lion | Evolution of spatially structured host-parasite interactions[END_REF], 2016[START_REF] Penczykowski | Understanding the ecology and evolution of host-parasite interactions across scales[END_REF]. Typically, however, in the host-parasite literature dispersal is considered as a fixed parameter in theoretical models or set by the experimenter in empirical studies (e.g., [START_REF] Morgan | The effect of migration on local adaptation in a coevolving host-parasite system[END_REF]. Therefore, it is largely unknown how dispersal itself may evolve in response to interactions with parasites, and how this might feed back on to epidemiological and coevolutionary processes. Dispersal has long been recognised as an important driver of metapopulation dynamics [START_REF] Hanski | Metapopulation Ecology[END_REF][START_REF] Clobert | Dispersal Ecology and Evolution[END_REF], and over recent years there has been an increasing interest in the evolution of dispersal [START_REF] Bowler | Causes and consequences of animal dispersal strategies: relating individual behaviour to spatial dynamics[END_REF][START_REF] Ronce | How does it feel to be like a rolling stone? Ten questions about dispersal evolution[END_REF][START_REF] Saastamoinen | Genetics of dispersal[END_REF]. It is now well established that dispersal has a genetic basis [START_REF] Saastamoinen | Genetics of dispersal[END_REF]). Thus, dispersal may readily evolve in response to spatial variation in fitness expectations [START_REF] Bowler | Causes and consequences of animal dispersal strategies: relating individual behaviour to spatial dynamics[END_REF], resulting from habitat variability [START_REF] Holt | Chaotic population dynamics favors the evolution of dispersal[END_REF], kin competition [START_REF] Hamilton | Dispersal in stable habitats[END_REF], or inbreeding depression [START_REF] Bengtsson | Avoiding Inbreeding -At What Cost?[END_REF][START_REF] Pusey | Inbreeding avoidance in animals[END_REF][START_REF] Perrin | Local competition, inbreeding, and the evolution of sex-biased dispersal[END_REF]. Only recently, theoretical work has begun to address dispersal evolution among interacting species in metacommunities (e.g., [START_REF] Poethke | Predator-Induced Dispersal and the Evolution of Conditional Dispersal in Correlated Environments[END_REF], Chaianunporn and Hovestadt 2012a,b, Amarasekare 2016). In the case of interactions with natural enemies, such as predators, herbivores or parasites, one may expect evolution of both constitutive and plastic dispersal behaviour. For example, constitutive longdistance dispersal in plants may evolve as a means to escape from specialised, locally dispersing pests and herbivores [START_REF] Muller-Landau | Theoretical perspectives on evolution of long-distance dispersal and the example of specialized pests[END_REF]. More recently, Chaianunporn and Hovestadt (2012b) showed that, if predation increases prey population extinction risk, higher constitutive dispersal rates may evolve to compensate this risk. In addition to such constitutive evolutionary changes, host dispersal behaviour may have evolved to be plastic in two main ways [START_REF] Boulinier | Dispersal and parasitism[END_REF]; see also [START_REF] Clobert | Informed dispersal, heterogeneity in animal dispersal syndromes and the dynamics of spatially structured populations[END_REF]): Dispersal may depend on the internal state of the host (state-dependence), that is, being infected or not, for example. In addition, dispersal may depend on the external context (context-dependence) which includes host densities or parasite prevalences, to name but two examples. Examples for the first case, that is, dispersal being infection state-dependent, have been reported for various organisms, from protists to vertebrates [START_REF] Sorci | Maternal parasite load increases sprint speed and philopatry in female offspring of the common lizard[END_REF][START_REF] Heeb | Ectoparasite infestation and sex-biased local recruitment of hosts[END_REF][START_REF] Suhonen | Activation of the immune system promotes insect dispersal in the wild[END_REF][START_REF] Fellous | Parasitic infection reduces dispersal of ciliate host[END_REF][START_REF] Brown | The impact of lungworm parasites on rates of dispersal of their anuran host, the invasive cane toad[END_REF]. However, the adaptive significance of these changes is not always clear. Reduced dispersal of infected hosts (e.g., [START_REF] Fellous | Parasitic infection reduces dispersal of ciliate host[END_REF] may simply be a physiological by-product of exploitation of host resources by the parasite. Conversely, increased dispersal may reflect manipulation by the parasite to increase parasite dispersal (Lion et al. 2006, Kamo and[START_REF] Kamo | The evolution of parasite dispersal, transmission, and virulence in spatial host populations[END_REF], rather than evolution of the host. Evolution of host dispersal plasticity has been studied theoretically by [START_REF] Iritani | Parasite infection drives the evolution of state-dependent dispersal of the host[END_REF] and [START_REF] Iritani | How parasite-mediated costs drive the evolution of disease state-dependent dispersal[END_REF]. They find that 'infectionbiased' dispersal, that is, increased dispersal when infected, is favoured with increasing probability of recovery during dispersal, with decreasing dispersal costs, and with higher levels of parasite virulence. Kin selection is the driving force underlying these predictions because the dispersal of infected individuals reduces the negative effects of kin competition on philopatric susceptibles, similar to [START_REF] Hamilton | Dispersal in stable habitats[END_REF]. In the latter case of dispersal plasticity, that is, context-dependent dispersal of hosts, dispersal occurs in response to the presence of natural enemies in a population, and the selective advantage is the possibility of escaping from mortality, just like in the above models of constitutive dispersal evolution. Examples of context-dependent dispersal are known from plants and animals where the presence of herbivores or predators induces an increased production of long-distance seeds (e.g., de la Peña and Bonte 2014) or a higher proportion of winged morphs in insects (e.g., [START_REF] Weisser | Predator-induced morphological shift in the pea aphid[END_REF]. This type of informed dispersal requires the capacity to 'sense' the presence of the enemy, either through direct contact when leaves are eaten by herbivores, for example, or indirect cues which can be mediated by chemicals released by predators or parasites or by damaged and infected hosts, for instance. [START_REF] Fronhofer | Bottom-up and top-down control of dispersal across major organismal groups[END_REF] showed that such chemical predator-related cues can indeed induce increased dispersal in various invertebrate and vertebrate organisms.

Accepted Ar ticle

'This article is protected by copyright. All rights reserved.' Despite the above discussed empirical examples, we are not aware of theoretical treatments of the evolution of context-dependent dispersal in host-parasite systems. Results of predator-prey models indicate that conditions for escape dispersal to evolve may be restrictive [START_REF] Poethke | Predator-Induced Dispersal and the Evolution of Conditional Dispersal in Correlated Environments[END_REF][START_REF] Hammill | Conspecific density modulates the effect of predation on dispersal rates[END_REF]. Namely, emigration from a currently infected population may not have a selective advantage if the probability of encountering the parasite is similar in other populations, or if encounter probabilities are unpredictable altogether. Moreover, if infected individuals take the same dispersal decisions as uninfected ones, escape benefits may be further limited. Using a spatially explicit, individual-based host-parasite model of a sterilizing parasite, we investigate conditions for the evolution of infection state-dependent dispersal and prevalence-dependent dispersal, where individuals can 'measure' local infection prevalence, that is, the relative number of infected hosts. Such a capacity may be a reasonable assumption, if, for example, prevalence is positively correlated with the quantity of a chemical parasite signal present in the environment, or if other signals, such as visual cues can be used to identify infection. Evolutionarily stable levels of plasticity are emergent properties, resulting from epidemiological and evolutionary dynamics in the host-parasite metapopulation. Importantly, our modelling framework takes kin competition into account by default [START_REF] Poethke | The relative contribution of individual and kin selection to the evolution of density-dependent dispersal rates[END_REF]. We find that infection-biased dispersal evolves over a wide range of conditions, mainly driven by kin selection, as previously described by [START_REF] Iritani | Parasite infection drives the evolution of state-dependent dispersal of the host[END_REF]. Going beyond previous theory, we show that the evolution of prevalence-dependent dispersal requires high levels of spatio-temporal autocorrelation in local prevalences, which arise for instance under high parasite virulence in our model where we observe negative prevalence-dependent dispersal. Under these conditions, prevalence-dependent dispersal feeds back on ecological dynamics and stabilizes the host-parasite system by reducing population fluctuations and extinction risk. Especially at low virulences, prevalence-dependent dispersal can also be positive if infected individuals pay increased dispersal costs.

The host-parasite model

Model overview

In this individual-based SI model, infected host individuals ( ) carry the parasite, and those that are not infected are susceptible ( ). For simplicity, there is no variation in resistance, and recovery is not possible. The parasite is sterilizing which implies that virulence reduces fecundity. The individuals disperse, reproduce and interact in a grid-like landscape ( patches) where each patch is connected to its eight nearest neighbours. In order to avoid edge effects, edges are wrapped around so that the landscape forms a torus. There are no differences in patch quality, and each patch initially has individuals, half of which are infected implying an initial infection prevalence of . The host has discrete, nonoverlapping generations and reproduces sexually. Dispersal has a genetic basis [START_REF] Saastamoinen | Genetics of dispersal[END_REF], is natal which implies that individuals may disperse before reproducing, and dispersal is costly [START_REF] Bonte | Costs of dispersal[END_REF]. Parasite transmission is horizontal and parasite dispersal only occurs together with infected hosts. Mating, density regulation, transmission and dispersal decisions are all governed by local rules. Once per generation, patches may suffer form an external local patch extinction. Since we model host-parasite dynamics explicitly, prevalences are not fixed a priori in contrast to [START_REF] Iritani | Parasite infection drives the evolution of state-dependent dispersal of the host[END_REF]. We run simulations for a total of 10,000 time steps which is enough to reach (quasi) equilibrium (Fig. S8). Analyses are performed over 10 replicate simulation runs. We consider two types of dispersal plasticity: First, host dispersal may depend on the infection state (state-or condition-dependence; see [START_REF] Iritani | Parasite infection drives the evolution of state-dependent dispersal of the host[END_REF], that is, the dispersal probability of an individual changes according to whether it is infected or susceptible. Second, the dispersal probability may change according to the local infection prevalence (context-dependence). We therefore assume that the host individuals can estimate local prevalence, as the fraction of infected individuals in their patch without any error (see [START_REF] Bocedi | Uncertainty and the Role of Information Acquisition in the Evolution of Context-Dependent Emigration[END_REF], Poethke et al. 2016a,b for a treatment of uncertainty during information acquisition). The two types of dispersal plasticity can evolve independently.

Reproduction and inheritance

The local host population dynamics are modelled explicitly using a discrete-time population growth model, with density-regulation according to Beverton1957

(1)

Accepted Ar ticle

'This article is protected by copyright. All rights reserved.' where is the fecundity, and is the intraspecific competition coefficient (see Tab. 1 for an overview of parameters and tested values). The equilibrium density in the Beverton-Holt model is given by ^ . After dispersal, the individuals reproduce. Females choose a mate at random from their local patch and produce a number of offspring drawn from a Poisson distribution with a mean which models demographic stochasticity. The factor of two appears in order to keep interpretable at the population level despite half of the population being males which do not produce offspring themselves. In case of a parasite infection will be reduced as described below. Each offspring has an equal chance of being either male or female. Since we assume non-overlapping generations, all adults die after reproduction. An individual's dependence of dispersal probability on the local prevalence in it's natal patch (

) and on the individual's infection state (

) is genetically encoded. For convenience, we assume a cubic relationship between dispersal probability and local infection prevalence (technically the four coefficients correspond to four different loci; Eq. 2) that may differ depending on whether individuals are susceptible or infected (state-dependence). This allows the separate evolution of reaction norms for the infected and the susceptible status, and thus represents a simple way to model the simultaneous evolution of infectiondependent and prevalence-dependent evolution.

(2)

The choice of a cubic reaction norm represents a compromise between a function-valued trait approach that would not make any assumption on the shape of the dispersal reaction norm [START_REF] Dieckmann | Surprising evolutionary predictions from enhanced ecological realism[END_REF]Metz 2006, Fronhofer et al. 2015) and assuming a functional relationship with fewer parameters. The functionvalued trait approach is not feasible here, as preliminary tests showed that the genetic algorithm was not able to optimize efficiently the reaction norm due to the large number of parameters. Each of the eight modelled loci (four per reaction norm per state) has two alleles, which are inherited independently (that is, we assume no linkage), one from each parent. During inheritance, mutations may occur which is modelled by adding a random number drawn from a normal distribution to the parental allelic value with a probability , which is the mutation rate. Simulations begin with standing genetic variation present in the population, that is, the coefficients in Eq. 2 assigned to each individual are drawn from a uniform distribution between 0 and 1. If there is very little variation in local prevalence, for instance, due to the epidemiological dynamics that result from the specific parameter settings of a given simulation, selection on dispersal probability at prevalences that are usually not experienced by the individuals is extremely weak. The reaction norm at these prevalences can therefore not be adequately optimized by our genetic algorithm. Therefore, we only present results of evolutionary stable dispersal reaction norms for the prevalences observed in the respective simulation settings.

Infection costs and transmission dynamics

Parasites usually utilise host resources for their own growth and reproduction, making their presence costly for the host. In our model, this is captured by the virulence ( ) which is a reduction in the mean number of offspring of the host, if infected:

(3) with if and if . and are the maternal and paternal infection states, respectively. The realised number of offspring is drawn from a Poisson distribution with mean . Of course, parasites can also increase mortality. Due to our assumption of discrete, non-overlapping generations we follow the implementation of virulence as a reduction in fecundity as proposed by [START_REF] Chaianunporn | The role of mobility for the emergence of diversity in victimexploiter systems[END_REF]. One of the defining differences between predator-prey and host-parasite models is that infected hosts may still reproduce. As we assume that all adults die after reproduction due to non-overlapping generations, we would have to apply parasite-induced mortality before reproduction which would effectively make our model a predator-prey model (for context-dependent dispersal in predator-prey systems see, e.g., [START_REF] Poethke | Predator-Induced Dispersal and the Evolution of Conditional Dispersal in Correlated Environments[END_REF]. Our model and its results are therefore most likely applicable to 'sterilizing parasites' in the broad sense that are known from plants and invertebrates, for example (for a detailed discussion see O'Keefe and Antonovics 2002, [START_REF] Best | Resistance is futile but tolerance can explain why parasites do not always castrate their hosts[END_REF]. All host individuals are born uninfected and therefore susceptible to infection. The parasite is transmitted after the death of the parental generation. The infection probability for the offspring generation depends on the density of the infected individuals in the parental generation (

), and on the density of all the individuals in the offspring generation (

). As in Chaianunporn and Hovestadt (2012b), this is rep-

Accepted Ar ticle

'This article is protected by copyright. All rights reserved.' resented using a Nicholson-Bailey type model [START_REF] Nicholson | The Balance of Animal Populations. -Part I[END_REF] with a Holling Type II (Hol-ling1959a, b) functional response.

( ) (4)
Here, is the searching efficiency and is the transmission rate. While originally developed for hostparasitoid systems, searching efficiency in our host-parasite model can be seen as the likelihood with which a host individual becomes infected if it encounters a parasite propagule. The transmission rate then governs how many parasite propagules one infected host will release upon its death.

Dispersal

Dispersal is natal and individuals disperse with a certain probability , based on their infection state (

) and local prevalence ( ) of the disease in the natal patch, calculated according to Eq. 2, where each coefficient is the mean of both alleles at a given locus. The target patch is randomly chosen among the eight nearest neighbouring patches. As dispersal is costly [START_REF] Bonte | Costs of dispersal[END_REF], individuals experience mortality during dispersal at a rate . Since virulence impacts reproduction as detailed above, we do not assume that susceptible and infected individuals pay differential costs of dispersal in contrast to [START_REF] Iritani | Parasite infection drives the evolution of state-dependent dispersal of the host[END_REF]. We nevertheless explore the consequences of relaxing this assumption in the Supplementary Material (see Figs. S9-S11).

Analysing the effect of kin selection

In an individual-based metapopulation model, kin structure arises by default [START_REF] Poethke | The relative contribution of individual and kin selection to the evolution of density-dependent dispersal rates[END_REF]) and is known to influence the evolution of dispersal [START_REF] Hamilton | Dispersal in stable habitats[END_REF]. To assess the effect of kin selection on the evolution of prevalence-and infection state-dependent dispersal we performed additional simulations using the algorithm described in [START_REF] Poethke | The relative contribution of individual and kin selection to the evolution of density-dependent dispersal rates[END_REF], which removes all kin structure from the metapopulation. Briefly, before dispersal and after the transmission of the parasite, all the individuals are shuffled and redistributed across the metapopulation, while maintaining the pre-shuffling density, prevalence, and sex ratio in each patch. This allows all-else-being-equal comparisons of results with and without kin structure. However, in our model, we may have two kin selection-relevant effects: First, the effect of kin competition described above, which is well-known to be relevant for dispersal evolution. Second, kin structure can also lead to the infection of kin, which may additionally affect inclusive fitness. Importantly, breaking kin structure using the approach proposed by [START_REF] Poethke | The relative contribution of individual and kin selection to the evolution of density-dependent dispersal rates[END_REF] will prevent both effects. In order to better tease these apart, we analyse an additional scenario where, instead of shuffling individuals, we randomly shuffle local prevalences before infection while keeping all other properties unchanged. This will break the association between prevalence of the parental generation and of the offspring generation locally and statistically prevent kin infection while conserving kin competition.

Alternative scenario: density-dependent dispersal

In general, infection prevalence may be correlated with population density. For instance, high prevalence could imply low host densities due to virulence, and vice versa. Therefore, instead of using prevalence as a source of information to base dispersal decisions on, hosts may evolve to use host population density as information, if population density is less costly to assess. To explore this possibility, we ran additional simulations replacing prevalence-dependent reaction norms by density-dependent reaction norms of dispersal. The dispersal rate was then a function of the local density (

) and infection state ( ). [START_REF] Poethke | Evolution of density-and patch-size-dependent dispersal rates[END_REF] have shown theoretically that in discrete-time models, such as ours, density-dependent dispersal must follow a threshold function of the form

{

(5)

In these simulations the threshold density corresponding to each infection state ( ) can evolve, instead of the polynomial coefficients in Eq. 2. The individuals do not disperse if the local density is below the threshold, and the dispersal probability increases in a saturating way thereafter.

Accepted Ar ticle

'This article is protected by copyright. All rights reserved.'

Results and Discussion

Both infection state-and prevalence-dependent dispersal evolved in our simulations (Fig. 1 A-C). Among the different demographic and epidemiological parameters in our model (see Tab. 1 for an overview), parasite virulence was a main factor impacting the evolution of dispersal plasticity qualitatively (see Supplementary Material Fig. S1-S6). For illustration, we have chosen a representative scenario, showing evolutionarily stable levels of dispersal plasticity for low, intermediate and high virulence (

, respectively), while other parameters are kept constant (Fig. 1 A-C).

Virulence modulates evolution of dispersal plasticity.

We find a transition from mainly infection state-dependent dispersal at low to intermediate levels of virulence ( ; Fig. 1 A,B) to clearly infection state-and prevalence-dependent dispersal at high virulence ( ; Fig. 1 C). This pattern is not very sensitive to other model parameters, and holds across different values of transmission rate ( ), parasite searching efficiency ( ), host fecundity ( ) and dispersal mortality ( ; see also Supplementary Material Fig. S1-S6).

In our example (Fig. 1 A-C) and more generally across a wide range of explored parameter values, infection state-dependent dispersal always evolves to become -biased, that is, dispersal is induced when individuals become infected. In fact, uninfected individuals practically never disperse, whereas infected dispersal probability is consistently high. The -bias becomes maximal at high virulence, where infection leads to unconditional dispersal (Fig. 1 C). This finding of an -bias confirms the theoretical predictions made by [START_REF] Iritani | Parasite infection drives the evolution of state-dependent dispersal of the host[END_REF]. We recover these results although, unlike their model, we do not assume a difference in the cost of dispersal for susceptible and infected individuals, nor do we allow for recovery of the infected individuals during dispersal. In addition, [START_REF] Iritani | Parasite infection drives the evolution of state-dependent dispersal of the host[END_REF] implement virulence on competitive ability, which hints that our results may be more widely applicable than to the sterilizing parasites we model. Furthermore, [START_REF] Iritani | Parasite infection drives the evolution of state-dependent dispersal of the host[END_REF] held prevalences constant across space and time, thus precluding the possibility for the evolution of prevalence-dependent dispersal. By contrast, we implemented the possibility for prevalence dependency of dispersal to evolve. Moreover, infection prevalence was allowed to vary freely from patch to patch, as a consequence of the interaction between dispersal, transmission dynamics (determined by and ) and virulence. It is important to note that the patterns shown in Fig. 1 may change when infected hosts pay higher costs of dispersal than susceptible hosts (Fig. S9). When the parasite has a very weak effect (low virulence) increasing dispersal costs for infected hosts trivially reduces -bias in dispersal. Reduced dispersal by infected hosts leads to more dispersal in susceptible hosts which becomes positively prevalencedependent. With increasing virulence the basic pattern described in Fig. 1 becomes more robust. -bias and negative prevalence-dependent dispersal for susceptible hosts at high virulences vanishes only if we assume that infected hosts pay a 3-fold increased dispersal cost in comparison to susceptible hosts (Fig. S9). In conclusion, a clear pattern of prevalence-dependent dispersal only evolves at high virulence and only for the susceptible state. Uninfected individuals show a steep negative reaction norm, meaning that dispersal is high when prevalence is locally low, but rapidly decreases to zero, as soon as prevalence reaches a threshold (Fig. 1 C). Infected individuals, by contrast, only show moderate trends of plasticity at low and intermediate virulence (Fig. 1 A,B). Below, we will discuss two mechanisms driving the observed evolution of dispersal plasticity.

Kin selection drives infection state-dependent dispersal.

Kin competition can be a major driver of dispersal evolution [START_REF] Bowler | Causes and consequences of animal dispersal strategies: relating individual behaviour to spatial dynamics[END_REF], even if dispersal is very costly [START_REF] Poethke | The relative contribution of individual and kin selection to the evolution of density-dependent dispersal rates[END_REF] or in the absence of ecological conditions favouring dispersal [START_REF] Hamilton | Dispersal in stable habitats[END_REF]. In our model, kin selection provides an intuitive explanation of -biased dispersal: infected individuals leave a patch to reduce the effect of kin competition on their uninfected relatives. This is because infected individuals have a lower fitness expectation than susceptible individuals. Therefore, the cost of dispersing in terms of inclusive fitness is lower if infected individuals disperse at a higher rate and susceptible individuals do not disperse. Hence, the dispersal rate evolves in such a way that kin competition is reduced (similar to [START_REF] Hamilton | Dispersal in stable habitats[END_REF] by being infection state-dependent. Indeed, both [START_REF] Poethke | Predator-Induced Dispersal and the Evolution of Conditional Dispersal in Correlated Environments[END_REF] for host-parasitoid and [START_REF] Iritani | Parasite infection drives the evolution of state-dependent dispersal of the host[END_REF] and [START_REF] Iritani | How parasite-mediated costs drive the evolution of disease state-dependent dispersal[END_REF] for hostparasite systems demonstrate that relatedness underlies the evolution of dispersal plasticity in the attacked organism. In addition to kin competition, a second mechanism may be relevant in our simulations:biased dispersal may also reduce the likelihood of kin infection. Existing theory has shown that even suicide may evolve as a strategy to reduce infection [START_REF] Débarre | Evolution of Host Life-History Traits in a Spatially Structured Host-Parasite System[END_REF].

Accepted Ar ticle

'This article is protected by copyright. All rights reserved.'

In individual-based metapopulation models like ours, kin structure arises by default [START_REF] Poethke | The relative contribution of individual and kin selection to the evolution of density-dependent dispersal rates[END_REF]). Its effect on evolutionary outcomes can be tested by using additional simulations with a constantly reshuffled genetic population structure [START_REF] Poethke | The relative contribution of individual and kin selection to the evolution of density-dependent dispersal rates[END_REF]. We find that this elimination of kin structure prevents the evolution of a clear dispersal bias at low and intermediate virulences, and thus susceptible and infected individuals disperse with similar probabilities for a given prevalence and virulence (Fig. 1 G,H). At high virulence (Fig. 1 I) the loss of context-dependency due to removed kin interactions secondarily leads to the extinction of the parasite due to very low dispersal rates (for a more comprehensive overview see Fig. S7). Both, the reduction of kin competition and of kin infection risk play a role in our simulations (Fig. 1). Avoiding kin infection alone seems not to be enough to explain the results we find, as the basic patterns can to some degree be found in simulations that prevent kin infection but not kin competition (Fig. 1 D-F). While the relative effects are difficult to quantify due to changing context-dependencies (Figs. S9-S11) the importance of kin infection versus kin competition depends on virulence: If one considers overall emigration rates, the effect of kin competition reduction becomes stronger with increasing virulence.

At low virulence (

; values correspond to results shown in Fig. 1 A, D and G) overall emigration rates are approximately 0.35 in the standard scenario, preventing kin infection reduces emigration rates to 0.22 and preventing kin competition and infection yield emigration rates of 0.1. At intermediate virulences (

; values correspond to results shown in Fig. 1 B, E and H) the corresponding rates are 0.73, 0.57 and 0.11. At high virulences removing kin competition has such a strong effect in reducing dispersal rates that the parasite ultimately goes extinct (Fig. 1 I). This effect of virulence can be readily explained because increased virulence reduces local population sizes (Fig. S14; see also Fig. 3) which in turn increases relatedness and kin competition. Vice versa, kin infection avoidance becomes less important at low virulences because of the reduced cost of infection. Unlike [START_REF] Iritani | Parasite infection drives the evolution of state-dependent dispersal of the host[END_REF], where -biased dispersal is partly explained by the opportunity for recovery during dispersal, there is no additional direct advantage of dispersal for infected individuals in our model. Thus, given the extreme effect of excluding kin structure, a plausible explanation for our -bias is that the dispersal of infected individuals increases the inclusive fitness of the philopatric susceptible individuals by reducing the effect of kin competition and kin infection. A priori, this risk should increase with virulence, and consequently the evolutionarily stable -bias should also increase with virulence, which is indeed what we observe in Fig. 1. One surprising finding is that dispersal bias sometimes evolves at no virulence and a non-zero dispersal mortality, that is, when there is no difference in fecundity between susceptible and infected individuals. Without any external extinction risk ( ), -biased dispersal evolves for zero and low virulences (Fig. S2). While counter-intuitive at first, we suggest that this is linked to an effect of demographic stochasticity [START_REF] Travis | The evolution of dispersal in a metapopulation: a spatially explicit, individual-based model[END_REF] acting on susceptible individuals: For high enough searching efficiencies, low virulences lead to a high prevalence, which means that the densities of susceptible individuals are greatly reduced. If an individual is susceptible in this scenario, it is more likely to go extinct. This interpretation is supported by simulations with greatly increased equilibrium densities ( , which implies an equilibrium density of 3000 individuals) where ESS dispersal rates tend towards 0 for both states, infected and non-infected.

Context-dependent dispersal evolves if conditions are sufficiently predictable.

Context-dependent dispersal only evolves if there is sufficient spatio-temporal variability in the relevant local condition (Poethke et al. 2016b). More specifically, plasticity is strongly favoured if spatial variation is combined with temporal autocorrelation in local conditions, that is, if conditions are sufficiently predictable for individuals to take the appropriate dispersal decision. We inspected the spatio-temporal variation in infection prevalence, the 'local context' in our model, across virulences. At low virulence ( ), prevalence is constantly high and does not show variation between patches (Fig. 2 A), which implies that context-dependent dispersal does not evolve because the context (prevalence) is not variable enough. The range of prevalences broadens for intermediate virulence and is maximal for high virulence (Fig. 2 A), with prevalences ranging from 0 to 1 across the metapopulation. Moreover, only for high virulence, we find strong oscillations in prevalence within patches, leading to generally higher levels of temporal autocorrelation (positive or negative; Fig. 2 B-D). This combines maximum spatial variation with temporal predictability, both prerequisites for dispersal plasticity to evolve. In other words, 'knowing' that low current infection prevalences will predictably change to high levels, explains the evolution of increased dispersal of susceptible individuals at low prevalence (Fig. 1 C), because this decision can be expected to yield an inclusive fitness benefit. Conversely, the rule not to disperse at high local prevalences means that staying in the local patch and waiting for better conditions (i.e., lower prevalences) is more beneficial than leaving the patch.

Accepted Ar ticle

'This article is protected by copyright. All rights reserved.'

The same arguments explain the evolution of a positive prevalence-dependent reaction norm for susceptible hosts when dispersal costs are higher for infected hosts: The -bias is lost and dispersal is generally low for infected hosts which become largely philopatric. Infected hosts staying in their natal patch make prevalences correlated in time, therefore uninfected individuals should disperse more when prevalences are high in order to potentially escape from infection. We observe that when there is no external extinction risk (

) and low parasite searching efficiency ( ), even high virulences ( ) may not lead to the evolution of prevalence-dependent dispersal, even though the entire range of prevalences occur (Fig. S2). Under these conditions, we observe a positive temporal autocorrelation in prevalences within a patch (see Fig. S12). Therefore, susceptible individuals stay in their natal patch at low prevalences because it is likely to remain that way giving them an inclusive fitness advantage. The spatial variation in prevalences implies that if individuals disperse, they are more likely to reach a patch with a higher prevalence, causing a greater infection risk for their offspring which reduces their fitness expectation [START_REF] Hastings | Can spatial variation alone lead to selection for dispersal[END_REF]. As the ESS emigration rate depends on the availability of good habitat in the metapopulation [START_REF] Poethke | A metapopulation paradox: partial improvement of habitat may reduce metapopulation persistence[END_REF] this selects against dispersal at low prevalence.

Evolution of infection state-and density-dependent dispersal.

Our model uses infection prevalence for context-dependent dispersal. However this information may not always be easily obtained by individuals. Population density is another, potentially more accessible cue that can drive context-dependent dispersal [START_REF] Travis | Habitat persistence, habitat availability and the evolution of dispersal[END_REF][START_REF] Poethke | Evolution of density-and patch-size-dependent dispersal rates[END_REF][START_REF] Matthysen | Density-dependent dispersal in birds and mammals[END_REF][START_REF] Nowicki | Evidence for positive density-dependent emigration in butterfly metapopulations[END_REF]. Because host density is also likely to be sensitive to parasite infestation, we re-ran our model using the density-dependent dispersal function derived by [START_REF] Poethke | Evolution of density-and patch-size-dependent dispersal rates[END_REF] as a reaction norm keeping the possibility for the evolution of state-dependency of dispersal. As above, we find the evolution of a consistent -bias of dispersal (Fig. 3 A-C), showing that state-dependent dispersal evolves independently of the type of context-dependency. We also find density-dependent dispersal of susceptible individuals under high virulence (Fig. 3 C). But unlike for prevalence dependency (Fig. 1 C), the relationship is positive: dispersal increases with increasing population density (Fig. 3 C). These mirrored patterns can be explained by the fact that density and prevalence are negatively correlated (Fig. 3 D-F): host populations become more decimated at higher infection prevalence if virulence is sufficiently high. Obviously, just like for prevalence, density needs to fluctuate in a predictable way for plasticity to evolve, which only occurs in the high-virulence scenario (Fig. 3 F). This finding implies that hosts may use either density or prevalence as a cue for contextdependent dispersal, depending on which information is less costly to obtain or more reliably estimated. Even though density and prevalence are generally negatively correlated, the relationship may not always be that simple. Due to the dynamic nature of population regulation by the parasite, prevalence and density are not perfectly negatively correlated, such that high population density may at times be associated with high prevalence, but at others with low prevalence. We show here that distinct patterns of cross correlations may arise, that is, prevalence and density are correlated with a time lag (Fig. 3 D-F).

Eco-evolutionary dynamics

Context-dependent dispersal based on density or on prevalence cues may have different demographic consequences at the metapopulation level. For instance, at ESS, patch extinction frequency, population density and prevalence across time can be impacted differentially (Fig. 4). Under the high-virulence scenario, where we observe the evolution of prevalence-dependent dispersal (Fig. 1 C), prevalence-dependency feeds back on demography and leads to a lower overall patch extinction rate than density-dependent dispersal, and to a higher frequency of intermediate prevalences at equilibrium (Fig. 4 and Figs. S13-S14). Therefore, while not reducing overall mean levels of prevalence, infection state-and prevalence-dependent dispersal reduces the occurrence of extremes, both in terms of prevalences and densities (Fig. 4). Prevalence-dependent dispersal is therefore stabilizing and may even be seen as a 'spatial tolerance' mechanism which can be relevant when the evolution of resistance is not possible. Investigating the concurrent evolution of resistance and other relevant life-history traits is of course an important next step. In summary, the demographic conditions defined by the host-parasite interaction drive the evolution of context-dependent dispersal which then feeds back on these demographic conditions.

General discussion and perspectives

In conclusion, we show that -biased, state-dependent dispersal readily evolves in our model and is pervasive. The dispersal bias evolves due to kin selection which is consistent with previous studies (Iritani and Accepted Ar ticle 'This article is protected by copyright. All rights reserved.' [START_REF] Iritani | Parasite infection drives the evolution of state-dependent dispersal of the host[END_REF][START_REF] Iritani | How parasite-mediated costs drive the evolution of disease state-dependent dispersal[END_REF]. Interestingly, kin selection processes can act both via kin competition reduction and via kin infection avoidance. Evolution of -biased dispersal is also possible, but only under increased dispersal costs for infected hosts (Fig. S9), which may arise when infection weakens the host and reduces dispersal capacity [START_REF] Binning | Parasites and Host Performance: Incorporating Infection into Our Understanding of Animal Movement[END_REF]. Alternatively, [START_REF] Iritani | Parasite infection drives the evolution of state-dependent dispersal of the host[END_REF] show that -bias may evolve, for example, if virulence reduces competitive ability or if parasite release occurs during dispersal, specificities not included in our model for the sake of simplicity. In contrast to state-dependent dispersal, we find rather restricted conditions for prevalence-, that is, context-dependent dispersal to evolve. Prevalence-dependent dispersal requires particular life-history trait values, such as high virulence, and a sufficient degree of predictability of epidemiological dynamics to evolve. In addition, the observed evolutionary outcome, a negatively prevalence-dependent dispersal reaction norm if susceptible, seems counter-intuitive at first. One may have naively predicted that uninfected individuals should only leave a patch if local prevalence is too high. In predator-prey systems, this seems to be a general pattern as recently demonstrated by [START_REF] Fronhofer | Bottom-up and top-down control of dispersal across major organismal groups[END_REF]. However, even in predator-prey systems, theory shows that the slope of the reaction norm will depend on the spatial and temporal autocorrelation of predation risk [START_REF] Poethke | Predator-Induced Dispersal and the Evolution of Conditional Dispersal in Correlated Environments[END_REF] as in our present study. Importantly, in our eco-evolutionary model, host-parasite dynamics and spatio-temporal heterogeneity in prevalences are a result of transmission which depends on the density of infected and susceptible individuals as well as on virulence and the evolved dispersal strategies. The emergent spatial variation in prevalences, temporal autocorrelation and oscillation in density drive the evolution of prevalence-dependent dispersal. Vice versa, the evolutionarily stable dispersal strategy modifies ecological dynamics. Specifically, prevalencedependent dispersal has a stabilizing effect (Fig. 4; see also [START_REF] Fronhofer | Bottom-up and top-down control of dispersal across major organismal groups[END_REF] for similar findings in food chains). Finally, we show that context-dependency in dispersal may evolve even if information on prevalences is not available, too costly to acquire or too noisy. Due to cross-correlations between prevalence and density, individuals may use local population density as an alternative source of information for contextdependent dispersal in host-parasite systems. Whether hosts can 'sense' the presence (and amount) of infection in their environment (e.g., via chemical cues) or simply use density as a cue, may be testable in experimentally tractable systems [START_REF] Fellous | Quorum Sensing and Density-Dependent Dispersal in an Aquatic Model System[END_REF][START_REF] Fronhofer | Bottom-up and top-down control of dispersal across major organismal groups[END_REF]. Empirical studies often describe negative (e.g., [START_REF] Cameron | Effects of body size and parasite infection on the locomotory performance of juvenile toads, Bufo bufo[END_REF][START_REF] Bradley | Parasites hinder monarch butterfly flight: implications for disease spread inmigratory hosts[END_REF][START_REF] Fellous | Parasitic infection reduces dispersal of ciliate host[END_REF][START_REF] Binning | Parasites and Host Performance: Incorporating Infection into Our Understanding of Animal Movement[END_REF][START_REF] Nørgaard | Infection in patchy populations: Contrasting pathogen invasion success and dispersal at varying times since host colonization[END_REF]), but also neutral [START_REF] Brown | The impact of lungworm parasites on rates of dispersal of their anuran host, the invasive cane toad[END_REF] or even positive effects [START_REF] Suhonen | Activation of the immune system promotes insect dispersal in the wild[END_REF][START_REF] Martini | Infection of an insect vector with a bacterial plant pathogen increases its propensity for dispersal[END_REF] of infection on dispersal-related traits. However, the adaptive significance of these observations is far from being clear (but see [START_REF] Tschirren | Maternal Modulation of Natal Dispersal in a Passerine Bird: An Adaptive Strategy to Cope with Parasitism?[END_REF], and may in some cases even reflect parasite manipulation [START_REF] Binning | Parasites and Host Performance: Incorporating Infection into Our Understanding of Animal Movement[END_REF][START_REF] Lion | The evolution of parasite manipulation of host dispersal[END_REF]). Generally, still few studies exist on context-dependent dispersal in host-parasite systems (Boulinier et al. 2001, French and[START_REF] French | Density-dependent dispersal in host-parasitoid assemblages[END_REF]. Our model was not specifically tailored to a particular host-parasite system. Instead we considered a relatively simple and tractable model, which nevertheless captures important features of natural systems. We aimed at a first step towards a general framework for plasticity evolution, taking into account the action of very general evolutionary processes, such as kin selection. Importantly, we extended previous work by relaxing strong assumptions of constant global prevalence, no transmission parameters and fixed population dynamics (see e.g., [START_REF] Iritani | Parasite infection drives the evolution of state-dependent dispersal of the host[END_REF]. This allows for eco-evolutionary feedbacks [START_REF] Govaert | Eco-evolutionary feedbackstheoretical models and perspectives[END_REF] to occur in our model (Fig. 4). Of course, the model can be extended and applied to more complex host-parasite systems, e.g., by including vertical transmission or different modes of sexual reproduction. Future challenges will be the explicit modelling of parasite life-history and (co)evolution of dispersal plasticity together with other traits, such as resistance or virulence.

Acknowledgements:

The authors would like to thank two anonymous referees for very helpful and constructive comments on previous versions of the manuscript. This is publication ISEM-YYYY-XXX of the Institut des Sciences de l'Evolution -Montpellier.

Funding: This study was funded by the CNRS.

Code availability: The computer code will be published using a Zenodo DOI upon acceptance.

Author contributions: J.N.D., O.K. and E.A.F. conceived the study. J.N.D. and E.A.F. developed the models. J.N.D. analysed the models. J.N.D. and E.A.F wrote the manuscript and all authors commented on the draft.

Accepted Ar ticle

'This article is protected by copyright. All rights reserved.'

Accepted Ar ticle

'This article is protected by copyright. All rights reserved.' for A, D and G, for B, E and H, and for C, F and I. This graph shows the reaction norms of 100 randomly chosen individuals at the end of the simulation (points) weighted by the normalised frequency of prevalences across space and time, that is, darker shades correspond to higher frequency of prevalence and lighter to lower frequency across space and time. We chose this weighting because the ecological settings of the simulations dictate the prevalences that occur and therefore also whether selection can act on parts of the reaction norm. Our approach (weighting combined with total simulation time, see Fig. S8) therefore guarantees that the results we show are ESS and not the result of drift, for instance, which will dominate for sections of the reaction norm that correspond to prevalences which have never occurred during a given simulation. The lines represent the median dispersal probability corresponding to that prevalence and the shaded regions are the quartiles. Note that the parasite goes extinct in I. Constant model parameters: , , , , .

Accepted Ar ticle

'This article is protected by copyright. All rights reserved.'

Accepted Ar ticle

'This article is protected by copyright. All rights reserved.'

Accepted Ar ticle

'This article is protected by copyright. All rights reserved.'

The distribution (frequency) of number of extinct patches (due to external local patch extinctions and emergent from the host-parasite dynamics), host densities and median prevalences at a given time for prevalence-dependent dispersal (PDD; red), density-dependent dispersal (DDD; orange) and an ecological scenario (RD; purple), for the last 1000 time steps. The ecological scenario is shown as a comparison and assumes that the dispersal rate is constant but equal to the observed overall ESS dispersal rate in the prevalence-dependent case. from left to right for the focal scenario ( , , , and , ). As expected, the results for and are very similar between the dispersal scenarios as under these conditions context-dependent dispersal does not evolve. For , the main effect of prevalence-dependent dispersal is a stabilizing effect: prevalence-dependent dispersal reduces the number of extinctions as well as extremes in host density and in prevalence.

Accepted Ar ticle

'This article is protected by copyright. All rights reserved.'


Ksenia Kolobova

William Rendu

Alena Shalagina

Pavel Chistyakov

Vasilii Kovalev

Malvina Baumann

Anastasia Koliasnikova

Andrei Krivoshapkin

The application of geometric-morphometric shape analysis to Middle Paleolithic bone retouchers

Keywords:

  L'archive ouverte pluridisciplinaire

Introduction

The Paleolithic of the Altai Mountains in southern Siberia is now widely known, not only because of the recent paleogenetic studies [START_REF] Slon | The genome of the offspring of a neanderthal and a denisovan father[END_REF], but also due to the discovery of the earliest known ornamental and bone tool complex in Eurasia, dating to the transition from the terminal Middle to Initial Upper Paleolithic [START_REF] Krivoshapkin | Between Denisovans and Neanderthals: strashnaya cave in the Altai Mountains[END_REF][START_REF] Douka | Age estimates for hominin fossils and the onset of the upper palaeolithic at denisova cave[END_REF]. Middle Paleolithic bone tools, on the other hand, have not been documented until recently [START_REF] Kolobova | Bone retouchers in the middle paleolithic complexes of Chagyrskaya cave[END_REF][START_REF] Baumann | The middle and upper paleolithic bone tools from the south chamber of denisova cave[END_REF][START_REF] Kozlikin | Bone retouchers from the denisova cave: new evidence[END_REF]. Two waves of Neanderthal in-migration to the Altai region have been recently defined on the basis of genetic and archaeological evidence. The origin of the Chagyrskaya Cave Neanderthals has been traced back to Central-Eastern Europe and the Caucasus region. A second Neanderthal migration took place around 70-60 kyr. B.P., bringing to the Altai a specific lithic assemblage similar to the European Micoquian/Keilmessergruppen characterized by the presence of radial cores and bifacial reduction sequences. This assemblage is predominated by simple and convergent side-scrapers and retouched points (Kolobova et al., 2020a, b). Bifacial plano-convex scrapers and points including Keilmesser type are numerous (144 items). These lithic artifacts occur in association with numerous bone tools (780 items), most of which are classified as retouchers (680 items, 87,2%) [START_REF] Baumann | Neandertal bone industry at Chagyrskaya cave (Altai, Russia)[END_REF]. The identification of these bone tool complexes raises numerous questions related to the cognitive abilities of Central Asian Neanderthals and Denisovans as well as their subsistence strategies, which have not yet been investigated. In this study, we focus on bone artifact assemblage from Chagyrskaya Cave, which has yielded the largest known osseous complex in the region [START_REF] Kozlikin | Bone retouchers from the denisova cave: new evidence[END_REF]. Recent advances in archaeological, paleoanthropological and paleontological researches have demonstrated the great scientific value and efficiency of geometricmorphometric analyses based on high-quality 3D models [START_REF] Otarola-Castillo | Differentiating between cutting actions on bone using 3D geometric morphometricsand Bayesian analyses with implications to human evolution[END_REF][START_REF] Chen | A late middle Pleistocene denisovan mandible from the Tibetan plateau[END_REF][START_REF] Herzlinger | A new method for 3D geometric morphometric shape analysis: the case study of handaxe knapping skill[END_REF][START_REF] Herzlinger | AGMT3-D: a software for 3-D landmarks-based geometric morphometric shape analysis of archaeological artifacts[END_REF][START_REF] Shalagina | Application of 3D geometric-morphometric analysis to the study of Stone Age lithic artifacts[END_REF]. Because bone retouchers are at the interface of two crucial activities for the subsistence of the groups: the butchery and the lithic tool production, the application of this analytical approach should provide new insights for our understanding of past human behaviors. However, due to the high degree of morphological variability among bone retouchers resulting from a combination of factors (bone anatomy, breakage patterns, taphonomy, etc.), geometric-morphometric approaches were not applied to these specific components of the material culture. In addition, sample size constraints (a minimum of 30 specimens, according to limitation of PCA) have strictly limited the technique's effective application. Here, we test the efficacy of geometric-morphometric analyses for quantitatively characterizing the retoucher shape. This study is a first step in the application of this method to the entire Altai Middle Paleolithic record.

Chagyrskaya Cave

Chagyrskaya Cave (51о26′34.6′′ N; 83о09′18.0′′ E) is situated on the left bank of the Charysh River in the Tigirek Range of the Altai Mountains of southern Russia. The 2007-2019 excavation campaign yielded significant lithic (117,421 specimens) and bone tool (780 items) assemblages as well as faunal and Neanderthal remains. The stratigraphy is composed of several archaeological layers (Layers 5-7) attributed to the final Middle Paleolithic. Layer 6 yielded the material used for this specific retoucher analysis. A series of absolute dates attribute the Neanderthal occupation of the cave to a relatively short period during final MIS4 and the beginning of MIS3. Paleoenvironmental data suggest that a steppe or semi-desert steppe environment prevailed under dry continental climatic conditions in the Charysh Valley at this time [START_REF] Derevianko | Multidisciplinary Studies of Chagyrskaya Cave -Middle Paleolithic Site in Altai[END_REF]Kolobova et al., 2020a,b). Abundant evidence of Neanderthal activities was uncovered in Layers 5-6c2. Layer 6c2, the lowermost in the Upper Pleistocene sequence, is preserved essentially in situ. The faunal assemblages from Layers 6a-6c2 include numerous bison remains (up to 47.2%) resulting from Neanderthal hunting activities that focused mainly on young individuals and females. The human impact on the Layer 6c2 assemblage is very high (31% of bones bear cut-marks, traces of scraping, and notches), while carnivore impact on the material is limited (8%) (for a synthesis of the paleontological and zooarchaeological data, see [START_REF] Derevianko | Multidisciplinary Studies of Chagyrskaya Cave -Middle Paleolithic Site in Altai[END_REF]Kolobova et al., 2020b). The systematic search for anthropogenic modification of faunal remains from Chagyrskaya Cave led to the identification of a new Middle Paleolithic bone tool assemblage composed of 680 retouchers. These artifacts are associated with diverse intermediate tools, retouched tools and implements with rounded tips, the functions of which remain to be determined [START_REF] Kolobova | Research of the middle paleolithic complexes from Chagyrskaya cave in 2017[END_REF][START_REF] Kolobova | The middle paleolithic complexes of the Chagyrskaya cave: 2018 field results[END_REF]Baumann et al., 2020, this volume). Chagyrskaya lithic assemblages have been associated with the Sibiryachikha Altai Middle Paleolithic variant, the easternmost manifestation of the Micoquian techno-complex (Kolobova et al., 2020a,b). Assemblages from Chagyrskaya Cave (layers 6a-6c2, 2008 excavation campaign) are characterized by flake production from radial and orthogonal cores; most flakes have trapezoidal or triangular shapes with relatively large flat or partly faceted striking platforms. The numerous toolkits (up to 30% of the assemblage, not including chips and debris) are dominated by simple and convergent scrapers of various shapes (70% of toolkits) and retouched points (15%). Bifacial flakes comprise up to 7% of the Chagyrskaya Cave lithic assemblage; a small component, but culturally significant in the Altai region. Here, Neanderthals undertook the selection of high-quality raw material to produce highly modified tools including convergent scrapers and bifacial implements [START_REF] Derevianko | Lithic raw material exploitation in the Sibiryachikha facies, the middle paleolithic of Altai[END_REF]Kolobova et al., 2019a). Retouchers define the interface of two technical systems; those directed toward animal exploitation and lithic production [START_REF] Costamagno | Bone retouchers and site function in the Quina mousterian: the case of les pradelles (Marillac-Le-France, France)[END_REF]; that comprise a formidable proxy for tackling the question of past hominin cognitive capacities through the analysis of a complex, multi-dimensional "chaine-op´eratoire." Layer 6c2 lithic assemblage demonstrate the frequent use of soft hammers by Neanderthals to produce unifacial and bifacial tools. The morphological characteristics of the proximal parts of flakes and core preparation spalls include numerous pronounced lips combined with diffused or absent bulbs in approximately 30% of the specimens examined. The pronounced bulbs and absent lips observed on the remaining flakes attest of the use of hard hammer percussion. Within the flake assemblage, only bifacial thinning flakes exhibit 100% of soft hammer reduction characteristics (Kolobova et al., 2019a). All these facts support a direct correlation between the numerous bone retouchers and the bifacial tools in Chagyrskaya Cave lithic assemblages. Prior to conduct our analysis, we have tested the efficiency of bone retouchers for the specific production of Chagyrskaya bifacial tools. According to our experimental model, the following technological sequences is observed: shaping of the flat surface, shaping of the convex surface and edge shaping. Our experiments attest of the greater efficiency of bone retouchers in comparison with hard and soft hammers for the final stages of the biface production. Moreover it demonstrates the advantages of bone retouchers in lithic tool retouching [START_REF] Shalagina | Application of 3D geometric-morphometric analysis to the study of Stone Age lithic artifacts[END_REF] for similar results see: [START_REF] Kolfschoten | Lower paleolithic bone tools from the 'spear horizon' at sch¨oningen (Germany)[END_REF]. The significant number of bone retouchers in the assemblage of Chagyrskaya Cave associated with the large number of unifacial tools strongly suggests that the first were also used for retouching the second.

3Materials and methods of analysis

Bone retouchers were sampled from a large assemblage (n = 680) derived from Chagyrskaya Cave cultural sub-layers 6a, 6b, 6c1 and 6c2. First of all, we excluded from the analysis specimens with major breakage in their functional areas (489 specimens, 72%) leaving 191 complete or almost complete retouchers, from which 40 specimens were randomly sampled from each layer for geometric-morphometric analyses. The entire sample was submitted to a thorough taphonomic analysis. Specimens were identified to the most precise taxonomic level possible and, when it was not possible to assign a specific attribution, ungulate size classes were used [START_REF] Brain | The Hunter or the Hunted. An Introduction to African Cave Taphonomy[END_REF]. Anatomical identification followed the "element, portion, segment" system of [START_REF] Gifford | A computer coding system for archaeological faunal remains[END_REF]. Bone surfaces of all specimens were first observed under low-oblique lighting using a 40X loupe for basic taphonomic and zooarchaeological observations. Evidence of weathering, root etching, anthropogenic and carnivore modification was systematically sought [START_REF] Olsen | Surface modification on bone: trampling versus butchery[END_REF][START_REF] Blumenschine | Blind tests of inter-analyst correspondence and accuracy in the identification of cut marks, percussion marks, and carnivore tooth marks on bone surfaces[END_REF][START_REF] Behrensmeyer | Taphonomic and ecologic information from bone weathering[END_REF]D'Errico and Villa, 1997;[START_REF] Pickering | Experimental patterns of hammerstone percussion damage on bones: implications for inferences of carcass processing by humans[END_REF]. Oxide coloration of bone cortical surfaces was also recorded and the proportion of preserved cortical surface was estimated per quartile [START_REF] Rendu | Hunting behavior and neanderthal adaptability in the late Pleistocene site of pech-de-l'Az´e I[END_REF][START_REF] Rendu | Subsistence strategy changes during the Middle to Upper Paleolithic transition reveals specific adaptations of Human Populations to their environment[END_REF]. When uncertain modification was detected, specimens were more thoroughly examined with a 20-80X light microscope. Overall, alteration was very rare; 90% of the bones exhibit complete or nearly complete preservation of their cortical surfaces. The principal taphonomic modifications result from light weathering or root etching. The assemblage of retouchers includes the following: The Chagyrskaya Cave bone tool assemblage includes a large number of complete or nearly complete retouchers (28%), 30% of the tools have post-depositional damages and 41% of the retouchers exhibiting marginal damages with dry fractures. The average length of complete retouchers -92.9 mm, width -37.2 mm, thickness -7.6 mm, weight -39.5 g. Our analysis focused on evaluating how this minor damages influenced retoucher shape and geometric-morphometric results. Complete retouchers without any intentional blank modifications reflect the variability of blanks and selection criteria. Complete retouchers exhibiting intentional edge modifications marked by spall negatives were analyzed to evaluate the selection of specific blank shapes by Paleolithic artisans. Our landmark-based geometricmorphometric shape analysis generated a quantitative description of variability within and between groups of specimens. This approach is based on the calculation of shape by measurement of the distances among landmarks and semi-landmarks located on the artifact, which were later processed in Cartesian space [START_REF] Herzlinger | A new method for 3D geometric morphometric shape analysis: the case study of handaxe knapping skill[END_REF]. We have employed Artifact GeoMorph Toolbox 3-D (AGMT3-D) to analyze retoucher shapes and followed the classical methodological approach of [START_REF] Herzlinger | AGMT3-D: a software for 3-D landmarks-based geometric morphometric shape analysis of archaeological artifacts[END_REF]. This analytical method has become widespread among archaeologists as an effective tool for evaluating the shape of Paleolithic artifacts [START_REF] Iovita | Comparing stone tool resharpening trajectories with the aid of elliptical fourier analysis. New perspectives on old stones: analytical approaches to paleolithic technologies[END_REF][START_REF] Iovita | Shape variation in aterian tanged tools and the origins of projectile technology: a morphometric perspective on stone tool function[END_REF][START_REF] Iovita | The handaxe reloaded: a morphometric reassessment of Acheulian and Middle Paleolithic handaxes[END_REF][START_REF] Archer | Diachronic change within the still bay at blombos cave, South Africa[END_REF][START_REF] Morales | The evolution and stability of stone tools: the effects of different mobility scenarios in tool reduction and shape features[END_REF]Herzlinger et al., 2017 etc.;[START_REF] Iovita | High handaxe symmetry at the beginning of the European Acheulian: the data from la Noira (France) in context[END_REF]. Generalized Procrustes Analysis (GPA) was used as a superimposition procedure, removing non-shape related variability, thus differences in landmark coordinates can be explained only by shape differences among retouchers. Thus, artifact size was not taken into account. The AGMT3-D package allows minimization of scaling by GPA, including minimizing the data dimension at each point by Least Squares criteria. Principal component analysis was employed for data dimensionality reduction and to provide a number of new axes (components) [START_REF] Lycett | A crossbeam co-ordinate caliper for the morphometric analysis of lithic nuclei: a description, test and empirical examples of application[END_REF][START_REF] Dryden | Holes and grooves: the contribution of microscopy and taphonomy to the problem of art origins[END_REF]. The analysis has included three steps: 3D scanning and modelling, positioning and measuring of landmarks and statistical processing of data. Digital 3-D models of retouchers were created using a RangeVision PRO 5M structured-light scanner. After scanning, models were processed using RangeVisionScanCenter and RangeVisionScanMerge software (Kolobova et al., 2019b). All models were then oriented in one perspective with the active retouching area in the northern or uppermost position. Positioning and measurement of landmarks were processed with AGMT3-D software. The landmarks were established within a dense grid of 30 × 30, resulting in 1800 landmarks recorded per artifact. According to our sample composition, we have defined three attributes corresponding to complete retouchers without blank modifications, complete retouchers with spall negatives and retouchers exhibiting only minor damage. In order to estimate the influence of the anatomical origin of bone blanks, we have defined six additional attributes according to the anatomical composition of the sample: femur, tibia, humerus, radius, vertebra and indeterminate long bone. We have performed statistical investigations with the relevant module in AGMT3-D and PAST-3 software. The PAST program was used for one-way parametric Multivariate Analysis of variance test (MANOVA) and a non-parametric PERMANOVA test was applied to the PCA scores analyses. The univariate non-parametric Wilcoxon Rank Sum Test was used to evaluate differences between samples [START_REF] Hammer | PAST: Paleontological statistics software package for education and data analysis[END_REF][START_REF] Herzlinger | AGMT3-D: a software for 3-D landmarks-based geometric morphometric shape analysis of archaeological artifacts[END_REF]. In order to describe retouchers exhibiting intentional spall negatives, we have performed a scar pattern/working steps analysis used to reconstruct tool-manufacturing sequences. Scar pattern analysis was based on the reconstruction of manufacturing and rejuvenation stages of lithic artifacts as well as on the chronological order of each stage in the manufacturing process [START_REF] Pastoors | Standardization and individuality in the production process of bifacial toolsleaf-shaped scrapers from the middle Paleolithic open air site Sare Kaya I[END_REF][START_REF] Kot | The earliest middle palaeolithic bifacial leafpoints in central and southern Europe: technological approach[END_REF][START_REF] Shalagina | Truncated-faceted pieces in the paleolithic of northern Asia[END_REF]Shalagina et al., 2019a,b).

Analytical results

Zooarchaeological characterization of the retoucher

Our sample is constituted of remains attributed to bison (NISP = 12), horse (NISP = 5) and non-identified large ungulates (NISP = 23) reflecting the other proportion already identified in the faunal spectrum. At the exception of one fragment of a transversal process of a bison vertebrae, all the blanks are made from long bones elements (Humerus: 5; Radius: 1; Femur = 3; Tibia = 15; unidentified long bone = 15). The preservation of the cortical surface is good (4 pieces have more than 50% of their cortical surface altered) and only the first stages of weathering are recorded on the retoucher sample. The only potential carnivore marks are pits on the bison vertebrae and a fragment of a humerus shaft, but because of some local exfoliation on their surface, it is not possible to exclude other possibilities. Fifty percent of the blanks exhibit evidence of scraping prior to their use as retouchers and 30% display notches linked to breakage of the bone for extracting marrow and grease.

Scar pattern analysis results

We have analyzed retouchers with negative scars mostly on cortical bone surfaces that could not have resulted from butchery or marrow consumption processes [START_REF] Grunwald | Analysis of fracture patterns from experimentally marrow-cracked frozen and thawed cattle bones[END_REF][START_REF] Morin | New criteria for the archaeological identification of bone grease processing[END_REF]. Recent experimental studies of marrow consumption at Chagyrskaya Cave reveal that during the flaking process it is virtually impossible to generate cortical negatives, in contrast with scars on the inner, medullary surfaces of bones [START_REF] Seletsky | Experimental Modelling of Bone Splitting Based on the Evidence of the Sibiryachikha Variant of the Altai Middle Paleolithic[END_REF]. The scar pattern analysis defines several types of spall negatives on bone surfaces, including those related to intentional blank modifications (Figs. 2,1234). The first cluster of negatives modified the initial bone shape (Figs.

2, 1, 2; negatives А). Some small negatives probably resulted from tool utilization (Figs. 2,1,negatives B,D). Several negatives used to sharpen retoucher edges are also observed (Figs. 2, 1, negatives Е; Fig. 2 and2, negatives С, Figs. 2, 3, negatives В; Figs. 2,4, negatives А, В). Since negative scars are located on various parts of bone blanks in a less regular way than in the case of stone artifacts, the complete reduction sequence is impossible to reconstruct, except when utilization traces overlap the main modification negatives (Figs. 2,1).

Geometric-morphometric analysis results

The PCA plot obtained with our geometric-morphometric approach present the variability of the retoucher shape (Fig. 4). The two first principal components account for 40.32% of the variability, which is a satisfactory outcome for geometric-morphometric approaches applied to paleoenvironmental reconstructions [START_REF] Kocovsky | The effect of sample size on the stability of principal components analysis of truss-based fish morphometrics[END_REF][START_REF] Shaukat | Impact of sample size on principal component analysis ordination of an environmental data set: effects on eigenstructure[END_REF]. The first three principal components account for more than 49.9% of the variability. These results demonstrate that all the retouchers analyzed have a relatively similar morphology. The variability of complete retouchers without blank modifications (193.27) is similar to the variability of retouchers bearing minor damage (193.2). To the contrary, retouchers bearing intentional spall negatives form a less variable group (166.15). All complete retouchers exhibiting intentional spall negatives are triangular in shape. The PERMANOVA multivariable test based on all PCA scores demonstrates no significant differences among the three retoucher groups (F: 0.6971; p: 0.8131). When all groups include several specimens, the sample composition allows us to use a parametric multivariable test. MANOVA, based on the three first PCA scores (F: 0.7477; p: 0.6133), confirmed the PERMANOVA results. We have performed paired tests to estimate the variability of the retoucher groups (Wilcoxon Rank Sum Test on within-group interpoint distances, on interpoint distances between-group means and on group centroid sizes). All the tests indicate an absence of difference between modified and unmodified retouchers on the one hand, and between unmodified retouchers and those with minor damages on the other. In order to understand the impact of initial shape of bone fragment on retoucher morphology, we have carried out a second geometric-morphometric analysis on the same sample. Among the 40 retoucher blanks analyzed, belonging to bisons and horses, 3 femora, 15 tibiae, 5 humeri, 1 radius, and 1 vertebra were identified, while 15 were assigned to the indeterminate long bone category. PCA plots generated from the geometric-morphometric analysis demonstrate the relative distribution of retouchers depending upon skeletal element attribution (Fig. 5). Because we analyzed the same sample, we observed the same variability coverage, thus we conclude that the nature of the blanks does not significantly influence the general variability of retouchers. At the same time, retouchers made from the same long bones demonstrate significant morphological differences (e. g., tibia or femur, Fig. 5). In the case of the small sample size groups (some groups are composed by only one specimen, such as a vertebra or radius), we were not able to use a parametric multivariable test. The PERMANOVA multivariable test based on all PCA scores demonstrates the similarity of all retouchers made on different blanks (F: 0.7782; p: 0.8732). If we exclude the two groups containing only one element (vertebra and radius), the MANOVA test based on the three first PCA scores (F: 0.5851; p: 0.8055) supports these results. Pair tests for estimating the variability of the retoucher groups (Wilcoxon Rank Sum Test on within-group interpoint distances, on interpoint distances between-group means and on group centroid sizes) also did not reveal any significant differences.

Discussion and conclusions

Appearing during the Lower Paleolithic, bone retouchers were frequently used during the Middle and Upper Paleolithic throughout the human-inhabited world. Eurasia has yielded the earliest evidence of bone blanks used as retouchers to fabricate lithic tools in Lower and Early Middle Paleolithic contexts: i.e., Caune de l'Arago [START_REF] Moigne | Utilisation des ossements de grands herbivores au Pal´eolithique inf´erieur[END_REF], La Micoque [START_REF] Langlois | Au sujet du Cheval de La Micoque (Dordogne) et des comportements humains de subsistance au Pl´eistoc`ene moyen dans le Nord-Est de l'Aquitaine[END_REF], Quesem Cave (Blasco et al., 2013b); Boxgrove [START_REF] Roberts | Boxgrove. A Middle Pleistocene hominid site at Eartham Quarry[END_REF], Terra Amata [START_REF] Moigne | Bone retouchers from Lower Palaeolithic sites Terra Amata, Orgnac 3, Cagny l'Epinette and Cueva del Angel[END_REF], and the Gran Dolina [START_REF] Rosell | Bone as a technological raw material at the gran Dolina site (sierra de Atapuerca, burgos, Spain)[END_REF]. In East Asia, the earliest bone retouchers are found at the Chinese Middle Paleolithic Lingjing site in Xuchang, Henan dated to 105-125 kyr BP (Doyon, Li, Hao, d'Errico, 2018). Numerous anatomical elements (e.g., upper limb bones, metapodials, skulls, mandibles, ribs) of large vertebrate taxa, mostly subjected to the butchery process, were used to produce various blanks. During the Upper Pleistocene, human societies may have assigned specific symbolic implications to the acquisition and use of unusual osseous tool blanks such as the Neanderthal skulls at La Quina, France [START_REF] Verna | The earliest evidence for the use of human bone as a tool[END_REF], Ursus spelaeus long bones at Scladina, Belgium [START_REF] Abrams | When Neanderthals used cave bear (Ursus spelaeus) remains: bone retouchers from unit 5 of Scladina Cave (Belgium)[END_REF] or, as in the case of some Aurignacian sites in southwestern France, large carnivore (Panthera spelaea) canines [START_REF] Castel | Sur Os et Sur Dents : les « retouchoirs » Aurignaciens de La Ferrassie (Savignac-de-Miremont, Dordogne)[END_REF]. Retouchers have generally been of interest to archaeologists only in the context as functional implements and are described as being produced on bone blanks, sometimes without evidence of carcass part selection [START_REF] Armand | Les retouchoirs en os d'Artenac (couche 6c): perspectives arch´eozoologiques, taphonomiques et exp´erimentales[END_REF]. They have often been considered informal, opportunistic tools made on blanks randomly obtained during the process of marrow extraction [START_REF] Chase | Tool-making and middle paleolithic behavior[END_REF][START_REF] Vincent | L'outillage osseux au Paleolithique moyen: une nouvelle approche[END_REF][START_REF] Armand | Les retouchoirs en os d'Artenac (couche 6c): perspectives arch´eozoologiques, taphonomiques et exp´erimentales[END_REF][START_REF] Patou-Mathis | Fiche g´en´erale[END_REF][START_REF] Karavani´c | The Middle Paleolithic Percussion or Pressure Flaking Tools? The comparison of experimental and archaeological material from Croatia[END_REF]. According to such hypotheses, there should be no morphological similarity with the exception of function-based size: weight ratios. Consequently, geometric-morphometric analyses have many advantages such as facilitating morphological classification without reference to the size. Our application of geometric-morphometric analyses to retouchers was successful. The most variable groups are composed of unmodified retouchers and retouchers exhibiting only minor damages. Some damaged retouchers fall outside the 95% ellipse of the complete retoucher variability, indicating that taphonomy can have a major impact by occluding artifact original forms (Fig. 4). Thus, in order to reduce general variability, it is necessary to exclude damaged retouchers from the analytical sample. Based on geometric-morphometric data, it is nearly impossible to separate complete retouchers from those exhibiting minor damages, hence taphonomic studies should always be conducted prior any further investigation on bone retouchers. Complete retouchers bearing intentional spall negatives (Fig. 3) constituted the least variable group. All of the modified retouchers are triangular or sub-triangular in plan form, which underscores the necessity of conducting research on specific shapes obtained due to human modification. The scar pattern analysis reveals that the modified portions of retouchers (Figs. 2,1,2) could significantly influence the geometric-morphometric results. The complete retouchers without blank modification constitute the most variable group. It is important to note that complete retouchers exhibiting intentional spall negatives fall within the middle range of retouchers without blank modifications 95% ellipse (Fig. 3). We assume that the blank modification was a significant technological stage in the retoucher shaping and the application of the approaches outlined in this study on a larger sample should confirm this hypothesis. Additionally, we have analyzed the influence of the anatomical origin of bone blanks. Our results (Fig. 5) do not demonstrate any significant effect of the anatomical origin of the bone retoucher blanks on their relative ordinal positions. Retouchers made from the same carcass elements often exhibit significantly variable shapes, suggesting that original bone morphology did not influence the morphological variability of the retouchers in our sample. This may indicate that while Neanderthals applied strict criteria to the selection of the retoucher blanks based on their physical properties (shape, weight, curvature, etc.), their specific anatomical origin was not taken into account. This may indicate the significant influence of intentional modification of blanks, whereby blanks of various shapes derived from different bones were intentionally modified to fall within certain parameters. In addition to define selectivity criteria of Middle Paleolithic retouchers [START_REF] Mallye | The mousterian bone retouchers of noisetier cave: experimentation and identification of marks[END_REF][START_REF] Auguste | Fiche ´eclats diaphysaires du Pal´eolithique moyen: biache-Saint-Vaast (Pas-de-Ca1ais) et Kulna (Moravie, R´epublique tch`eque)[END_REF][START_REF] Costamagno | Bone retouchers and site function in the Quina mousterian: the case of les pradelles (Marillac-Le-France, France)[END_REF], this could imply the existence of defined morphological patterns for retoucher blanks resulting from intentional selection and modification, which could change our understanding of Neanderthal subsistence strategies and cognitive abilities. There is a parallel between the predetermination of standardized shapes of lithic bifacial tools from Chagyrskaya Cave and the production of bone retouchers of a specific shape, accompanied by the absence of constraint on bone morphology. This implies that the selection and production of bone retouchers may have occurred very early in the butchery process, as [START_REF] Costamagno | Bone retouchers and site function in the Quina mousterian: the case of les pradelles (Marillac-Le-France, France)[END_REF] proposed and, consequently, would constitutes evidence of the very strong imbrication of carcass processing and the selection of bone tools for the production of lithic artifacts. This is additional evidence of the role of bone retouchers as an interface between lithic production and hunting activities, and highlights their great potential for addressing the question of predetermination in complex, multidimensional "chaine-op´eratoire" and, ultimately, discussions of the cognitive capacities of our hominin ancestors.

Declaration of competing interest

None.

  1. Complete retouchers without blank modifications (Figs. 1 and 14 specimens); 2. Complete retouchers with intentional blank modifications (Figs. 2 and 10 specimens); 3. Retouchers with minor damage (Figs. 3 and 16 specimens).

Acknowledgements

The excavation of Chagyrskaya Cave was funded by the Russian Foundation for Basic Research (project 18-09-00041). We are thankful to the joint Russian-French Scientific Project of the Russian Foundation for Basic Research (N• 19-59-22007) -Fondation Maison des sciences de l'Homme (FMSH) for support of the bone tool and geometric-morphometric analysis of this research. Paleontological analysis was supported by the Russian Foundation for Basic Research (N• 18-09-40070). Zooarchaeological analyses were funded by the CNRS, the Idex from Bordeaux University, and the French Ministry of Culture within the framework of the IRL Artemir. We gratefully acknowledge the anonymous reviewers and John Olsen for his valuable comments that deeply enriched the paper.