Léa Terray

Emmanuelle Stoetzel

David Nerini

Raphaël Cornette

How Can Phenotypic Evolution be Characterized Over Time and Through Environmental Changes?

Keywords: Phenotypic evolution, Paleoenvironments, Rodents, Teeth, Geometric morphometrics

Rapid phenotypic evolution is observed in response to rapid environmental changes. These phenotypic variations can occur at different scales, from the population to the community. We intended to characterize these multiscale phenotypic responses in rodents from the archaeological site El Harhoura 2 (Rabat, Morocco), dated from the Late Pleistocene to the Holocene, and relate them to paleoenvironmental variations. Upper and lower first molars were used as proxies for phenotype. Their shape was quantified using a landmark-free method. To account for both intra-and interspecific variations, we used morphogroups as phenotypic units. Those morpho-groups were identified using unsupervised clustering. Three shape indicators were computed: number of morphogroups, morpho-group disparity and changes in morpho-group mean shape over time. To evaluate this little-used approach on small mammals, we compared it to three widely used biodiversity indices: number of taxa, Shannon index, and Simpson index. Phenotypic evolution between and within species was highlighted. Morpho-groups seemed to be phenotypic response units representing ecological groups that transversed species. Variations in the morpho-group mean shapes were partly related to paleoenvironmental changes; however, variations in disparity were not. Thus, environmental changes deduced from fossil microvertebrate communities did not seem to be the main determinants of the characterized phenotypic variations.

Introduction

Climate change affects the living environments of organisms (vegetation, landscapes, etc.). These changes can be gradual or abrupt, and morphological variations at adaptative and/or phenotypical scales are often observed in response to these new environmental conditions [START_REF] Lande | Adaptation to an extraordinary environment by evolution of phenotypic plasticity and genetic assimilation[END_REF][START_REF] Hoffmann | Climate change and evolutionary adaptation[END_REF][START_REF] Clavel | Accelerated body size evolution during cold climatic periods in the Cenozoic[END_REF]. Rapid phenotypic evolution in response to changes in the environment is related to microevolution and phenotypic plasticity [START_REF] Holt | The microevolutionary consequences of climate change[END_REF][START_REF] Lande | Adaptation to an extraordinary environment by evolution of phenotypic plasticity and genetic assimilation[END_REF][START_REF] Hoffmann | Climate change and evolutionary adaptation[END_REF][START_REF] Boutin | Climate change and mammals: evolutionary versus plastic responses[END_REF][START_REF] Merila | Climate change, adaptation, and phenotypic plasticity: the problem and the evidence[END_REF], two processes that generate interindividual variations (e.g., [START_REF] De Jong | Phenotypic plasticity as a product of selection in a variable environment[END_REF][START_REF] Fusco | How many processes are responsible for phenotypic evolution?[END_REF]. However, the phenotype also results from macroevolution that generates variations at the interspecific level [START_REF] Erwin | Macroevolution is more than repeated rounds of microevolution[END_REF][START_REF] Hautmann | What is macroevolution?[END_REF]). Phenotypic evolution can thus refer differentially to intra-or interspecific evolution. The difference between those levels of variation is often not obvious because evolution is a continuous process, especially if there is an ongoing trend of divergence. Macroevolutionary changes are guided by microevolutionary changes occurring at each generation [START_REF] Hansen | Translating between microevolutionary process and macroevolutionary patterns: the correlation structure of interspecific data[END_REF], and the uncoupling between these two processes is still being debated [START_REF] Simpson | Tempo and Mode in Evolution[END_REF][START_REF] Erwin | Macroevolution is more than repeated rounds of microevolution[END_REF][START_REF] Arnold | The adaptive landscape as a conceptual bridge between micro-and macroevolution[END_REF][START_REF] Hendry | An introduction to microevolution: rate, pattern, process[END_REF][START_REF] Reznick | Darwin's bridge between microevolution and macroevolution[END_REF]. To characterize the phenotypic response to environmental changes, we chose to address phenotypic evolution as a composite process implying both intra-and interspecific evolution.

Most studies raising the issue of phenotypic evolution have used species as the phenotypic units to characterize phenotypic variation (e.g., [START_REF] Hendry | Human influences on rates of phenotypic change in wild animal populations[END_REF][START_REF] Bickford | Impacts of climate change on the amphibians and reptiles of Southeast Asia[END_REF][START_REF] Gardner | Declining body size: a third universal response to warming?[END_REF][START_REF] Reed | Interacting effects of phenotypic plasticity and evolution on population persistence in a changing climate: evolution, plasticity, and climate change[END_REF][START_REF] Boutin | Climate change and mammals: evolutionary versus plastic responses[END_REF][START_REF] Radchuk | Adaptive responses of animals to climate change are most likely insufficient[END_REF]. However, individuals belonging to the same species may display distinct phenotypic responses to environmental conditions, especially in a complex climatic context [START_REF] Sultan | Metapopulation structure favors plasticity over local adaptation[END_REF]. The concept of species itself, albeit discussed, defined as a morphologically or genotypically similar grouping of individuals, implies the existence of individuals with intermediate morphologies [START_REF] Mayr | Systematics and the origin of species[END_REF][START_REF] Simpson | The species concept[END_REF][START_REF] Mallet | A species definition for the Modern Synthesis[END_REF][START_REF] Wheeler | Species Concepts and Phylogenetic Theory: A Debate[END_REF][START_REF] Baker | Speciation in mammals and the genetic species concept[END_REF]. These intermediate morphotypes can display a variety of responses to environmental changes. When the addressed problematic goals are to identify patterns and/or characterize processes at a scale other than the species scale (as an intraspecific scale), this partition may not be appropriate. Therefore, in our case, that kind of approach is not the most suitable way to demonstrate phenotypic evolution on multiple scales. To detect phenotypic responses to environmental changes, we need to adopt a taxonomy-free approach.

In the context of our study, morphology is a potential candidate criterion to define those taxonomy-free units. Morphology as a phenotypic unit is particularly used in the field of study of foraminifera (e.g., [START_REF] Nagy | Environmental significance of foraminiferal morphogroups in Jurassic North Sea deltas[END_REF][START_REF] Tyszka | Response of Middle Jurassic benthic foraminiferal morphogroups to dysoxic/anoxic conditions in the Pieniny Klippen Basim Polish Carpathians[END_REF][START_REF] Alperin | Benthic foraminiferal morphogroups on the Argentine continental shelf[END_REF][START_REF] Alegret | End-Cretaceous marine mass extinction not caused by productivity collapse[END_REF][START_REF] Khare | Cluster analysis of benthic foraminiferal morphogroups from the western margin of India reflects its depth preference[END_REF]). This shape-based approach is based on setting morphological groups (MGs) rather than taxa as biodiversity units. This approach groups individuals solely on the basis of their morphology. In the context of foraminifera, [START_REF] Murray | Distribution and Ecology of Living Benthic Foraminiferids[END_REF][START_REF] Murray | Distribution and Ecology of Living Benthic Foraminiferids[END_REF] defined MGs as phylogenetically independent groups of forms with similar test morphology. Since a test's form is closely related to its environment, this MG approach allows one to infer climatic conditions [START_REF] Alperin | Benthic foraminiferal morphogroups on the Argentine continental shelf[END_REF]). Conversely, it can be a relevant way to characterize phenotypic responses to environmental changes when studied organisms are known to be good environmental indicators. The concept of MGs has been successfully extended to other study frameworks, such as soil biocrust in [START_REF] Read | Biocrust morphogroups provide an effective and rapid assessment tool for drylands[END_REF], where they concluded MGs were functional response groups to environmental disturbances. Moreover, the MG approach presents demonstrated advantages: 1) it allows rapid and ecologically informative surveys of large morphological datasets; 2) it enables comparisons of assemblages of different ages; 3) taxonomic identification is not needed; and 4) MGs are independent of systematic relationships [START_REF] Murray | Distribution and Ecology of Living Benthic Foraminiferids[END_REF][START_REF] Murray | Distribution and Ecology of Living Benthic Foraminiferids[END_REF][START_REF] Nagy | Environmental significance of foraminiferal morphogroups in Jurassic North Sea deltas[END_REF][START_REF] Read | Biocrust morphogroups provide an effective and rapid assessment tool for drylands[END_REF]. We extended Murray's definition and defined MGs as phylogenetic independent groupings of forms with similar morphology.

El Harhoura 2 (EH2), a coastal cave located in the Rabat-Temara region. This site is of particular interest for the study of short-term phenotypic evolution related to environmental variations [START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF][START_REF] Stoetzel | Adaptations and dispersals of anatomically modern humans in the changing environments of North Africa: the contribution of microvertebrates[END_REF][START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF][START_REF] Stoetzel | Taphonomy and palaeoecology of the late Pleistocene to middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF]Stoetzel et al. , 2012b[START_REF] Stoetzel | Systematics and evolution of the Meriones shawii/grandis complex (Rodentia, Gerbillinae) during the Late Quaternary in northwestern Africa: exploring the role of environmental and anthropogenic changes[END_REF]. This region is subjected to complex climatic influences from the Atlantic Ocean to the west, the arid Sahara to the south, and the Mediterranean region to the north. EH2 covers a time period from the Late Pleistocene to the Holocene, during which important climatic fluctuations occurred in the area (e.g., [START_REF] Drake | Ancient watercourses and biogeography of the Sahara explain the peopling of the desert[END_REF][START_REF] Drake | Palaeoclimate in the Saharan and Arabian deserts during the Middle Palaeolithic and the potential for hominin dispersals[END_REF][START_REF] Blome | The environmental context for the origins of modern human diversity: a synthesis of regional variability in African climate 150,000-30,000 years ago[END_REF][START_REF] Kageyama | Mid-Holocene and Last Glacial Maximum climate simulations with the IPSL model-part I: comparing IPSL_CM5A to IPSL_CM4[END_REF][START_REF] Scerri | The North African Middle Stone Age and its place in recent human evolution[END_REF][START_REF] Couvreur | Tectonics, climate and the diversification of the tropical African terrestrial flora and fauna[END_REF]. This variable climatic context resulted in a succession of relatively humid/arid and open/closed environments at EH2 [START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF][START_REF] Stoetzel | Taphonomy and palaeoecology of the late Pleistocene to middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF]Stoetzel et al. , 2012a, b), b). This cave revealed an exceptional richness in small terrestrial vertebrate remains accumulated through owl pellets and/or carnivore faeces, which have been intensively studied over the past decade (e.g., [START_REF] Michel | La macrofaune du Pleistocene superieur d'El Harhoura 2 (Temara, Maroc) : donnees preliminaires[END_REF][START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF][START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF][START_REF] Stoetzel | Taphonomy and palaeoecology of the late Pleistocene to middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF]Stoetzel et al. , 2012b[START_REF] Stoetzel | Systematics and evolution of the Meriones shawii/grandis complex (Rodentia, Gerbillinae) during the Late Quaternary in northwestern Africa: exploring the role of environmental and anthropogenic changes[END_REF][START_REF] Cornette | Shrews of the genus Crocidura from El Harhoura 2 (Temara, Morocco): the contribution of broken specimens to the understanding of Late Pleistocene-Holocene palaeoenvironments in North Africa[END_REF].

Among the small mammals of EH2, rodents are well represented in every layer. They are the most diverse group of living mammals [START_REF] Carleton | Order Rodentia[END_REF] and are known for being good paleoenvironmental indicators (e.g., [START_REF] Avery | Micromammals as paleoenvironmental indicators and an interpretation of the late Quaternary in the Southern Cape Province[END_REF][START_REF] Fernandez-Jalvo | Taphonomy and palaeoecology of Olduvai Bed-I (Pleistocene, Tanzania)[END_REF][START_REF] Mcguire | Geometric morphometrics of vole (Microtus californicus) dentition as a new paleoclimate proxy: shape change along geographic and climatic clines[END_REF][START_REF] Belmaker | Ecological change and the extinction of the Levantine Neanderthals: implications from a diachronic study of micromammals from Amud Cave, Israel[END_REF][START_REF] Arregoitia | Morphology captures diet and locomotor types in rodents[END_REF][START_REF] Royer | New bioclimatic models for the quaternary palaearctic based on insectivore and rodent communities[END_REF][START_REF] Lopez-Garcia | Rodents as indicators of the climatic conditions during the Middle Pleistocene in the southwestern Mediterranean region: implications for the environment in which hominins lived[END_REF]. Rodent morphology can thus reflect ecological similarities. Among most mammals, teeth display higher homoplasy than osteological parts [START_REF] Evans | High-level similarity of dentitions in carnivorans and rodents[END_REF][START_REF] Brocklehurst | Dental characters used in phylogenetic analyses of mammals show higher rates of evolution, but not reduced independence[END_REF]. Thus, rodent teeth may be a good proxy for estimating rodent phenotype. The shape-based grouping of fossil rodent incisors has been shown to reflect dietary ecology [START_REF] Paine | The ecomorphology of southern African rodent incisors: potential applications to the hominin fossil record (C Charles[END_REF], and molar crowns are strongly related to the environment through diet [START_REF] Wolf | Does weather shape rodents? Climate related changes in morphology of two heteromyid species[END_REF][START_REF] Mcguire | Geometric morphometrics of vole (Microtus californicus) dentition as a new paleoclimate proxy: shape change along geographic and climatic clines[END_REF][START_REF] Coillot | Correlated changes in occlusal pattern and diet in stem murinae during the onset of the radiation of old world rats and mice[END_REF][START_REF] Cano | Dietary ecology of murinae (Muridae, Rodentia): a geometric morphometric approach[END_REF][START_REF] Rodrigues | The great variety of dental structures and dynamics in rodents: new insights into their ecological diversity[END_REF][START_REF] Selig | Mammalian molar complexity follows simple, predictable patterns[END_REF].

We studied the phenotypic evolution of the rodents of the EH2 cave using an MG approach to answer the following question: does short-term phenotypic evolution occur at different scales, i.e., intra/interspecific, under changing environments? We used two complementary models: first lower molars (m1) and first upper molars (M1). In murine rodents, both display slightly different evolutionary rates of size and shape [START_REF] Renaud | Influence of biotic and abiotic environment on dental size and shape evolution in a Late Miocene lineage of murine rodents (Teruel Basin, Spain)[END_REF]. MGs were identified using a nonsubjective grouping method: unsupervised clustering [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]. For each dataset, three shape indicators were computed to evaluate phenotypic diversity: number of MGs, disparity over time and changes in mean shape over time. We expected MGs to not reflect species but to be representative of response groups towards environmental changes. We also expected to demonstrate both variations shared by several MGs under strong environmental changes (related to interspecific phenotypic evolution) and MG-specific variations when MGs display different phenotypic responses (related to intraspecific phenotypic evolution). The number of MGs and disparity should detect periods of decreasing/increasing/stable phenotypic diversity. Variations in mean shape should indicate shifts in selected shapes and the eventual convergence/divergence of MG shapes in constraining environments. However, the MG approach often underestimates diversity and can be costly in terms of information [START_REF] Read | Biocrust morphogroups provide an effective and rapid assessment tool for drylands[END_REF]. Then, to assess confidence in our approach, we compared its ability to characterize phenotypic evolution to widely used biodiversity indices in ecology and archaeology (number of species, Shannon index and Simpson index) (e.g., [START_REF] Avery | Micromammals as paleoenvironmental indicators and an interpretation of the late Quaternary in the Southern Cape Province[END_REF][START_REF] Geraads | The rodents from the late middle Pleistocene hominidbearing site of J'bel Irhoud, Morocco, and their chronological and paleoenvironmental implications[END_REF][START_REF] Meunier | Mise a jour de la liste des rongeurs d'Algerie, biogeographie et implications paleoecologiques[END_REF].

Material and Methods

El Harhoura 2 Cave

The stratigraphy of the cave is structured into 11 layers numbered from top to bottom (Fig. 1). Eight levels (L1 to L8) are ;well studied from their microfauna [START_REF] Stoetzel | Taphonomy and palaeoecology of the late Pleistocene to middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF][START_REF] Stoetzel | Mus in Morocco: a Quaternary sequence of intraspecific evolution[END_REF]2017) and well dated [START_REF] Jacobs | Single-grain OSL chronologies for Middle Palaeolithic deposits at El Mnasra and El Harhoura 2, Morocco: implications for Late Pleistocene human-environment interactions along the Atlantic Coast of northwest Africa[END_REF]Ben Arous et al. ;2020a, b;[START_REF] Marquer | The first use of olives in Africa around 100,000 years ago[END_REF] and were used in this study. A small mammal taxonomic study of the site reveals the presence of at least eight rodent species: Apodemus sylvaticus, Lemniscomys barbarus, Mus spretus, Dipodillus campestris, Meriones shawii, Meriones grandis, Jaculus cf. orientalis and Eliomys cf. munbyanus [START_REF] Stoetzel | Taphonomy and palaeoecology of the late Pleistocene to middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF][START_REF] Stoetzel | Mus in Morocco: a Quaternary sequence of intraspecific evolution[END_REF][START_REF] Stoetzel | Systematics and evolution of the Meriones shawii/grandis complex (Rodentia, Gerbillinae) during the Late Quaternary in northwestern Africa: exploring the role of environmental and anthropogenic changes[END_REF]. We were not able to include Jaculus cf. orientalis and Eliomys cf. munbyanus due to the lack of well-preserved teeth remains. Local paleoenvironmental information from EH2 is based on large and small vertebrates recovered from archaeological levels [START_REF] El Hajraoui | Prehistoire de La Region de Rabat-Temara[END_REF][START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF][START_REF] Stoetzel | Taphonomy and palaeoecology of the late Pleistocene to middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF]Stoetzel et al. , 2012b[START_REF] Stoetzel | Context of modern human occupations in North Africa: contribution of the Temara caves data[END_REF][START_REF] Campmas | Were Upper Pleistocene human/nonhuman predator occupations at the Temara caves (El Harhoura 2 and El Mnasra, Morocco) influenced by climate change?[END_REF][START_REF] Stoetzel | Adaptations and dispersals of anatomically modern humans in the changing environments of North Africa: the contribution of microvertebrates[END_REF]. They show successive relatively humid (L3, L4a, L6 and L8) and arid (L2, L5 and L7) phases. Late Pleistocene environments appeared more open and less humid than those today. Paleolandscapes are defined as open steppe or savanna-like land with patches of shrubland, woodland and water bodies. Water bodies cover greater areas during humid phases than during arid phases [START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF].

Data Collection

We used 1133 m1s and 990 M1s housed at the Musee de l'Homme, Paris, France and the Institut National des Sciences de l'Archeologie et du Patrimoine, Rabat, Morocco. The sample was composed of six of the eight species present at EH2: Meriones shawii, Meriones grandis (referred to in the following as the Meriones shawii/grandis complex because of their morphological proximity [START_REF] Stoetzel | Systematics and evolution of the Meriones shawii/grandis complex (Rodentia, Gerbillinae) during the Late Quaternary in northwestern Africa: exploring the role of environmental and anthropogenic changes[END_REF])), Mus spretus and Dipodillus campestris in abundance and occasionally Lemniscomys barbarus and Apodemus sylvaticus. The abundance of teeth per layer is indicated in Table 1.

Tooth photographs were taken using a Nikon digital camera D 5500 coupled with AF-S Micro NIKKOR 60 mm and macro extension tubes. Picture taking was standardized: teeth were always placed so that the occlusal surface was horizontal. We obtained two categories of pictures: in situ teeth (when teeth were set in the mandible or the maxilla) and individual teeth (when teeth were no longer encased in the bone). In the first case, masks of teeth were manually extracted using the software Gimp v 2.10.6 (The GIMP Development Team 2018). In the second case, mask extraction was semi-automatized using ImageJ 1.52j software [START_REF] Schneider | NIH Image to ImageJ: 25 years of image analysis[END_REF] through the use of the segmentation procedure of the MorphoLibJ plugin [START_REF] Legland | MorphoLibJ: integrated library and plugins for mathematical morphology with ImageJ[END_REF], which automatically detects objects on a picture. The difference between those two acquisition methods was tested and was not statistically significant (Online Resource 1). The data acquisition protocol is summarized in Fig. 2. Finally, we extracted outlines from the masks using the library Momocs [START_REF] Bonhomme | Momocs: outline Analysis Using R[END_REF] in the free software R (R Development Core Team 2018).

Geometric Morphometric Analyses

The tooth shape was approximated through 2D outlines of the teeth in occlusal view. Outlines were taken at the base of the molar crown, which is less affected by age. This is especially the case for Meriones, which has emihypsodont characteristics. As a consequence, the occlusal surface and pattern vary importantly with the age of the individual (as discussed in [START_REF] Stoetzel | Systematics and evolution of the Meriones shawii/grandis complex (Rodentia, Gerbillinae) during the Late Quaternary in northwestern Africa: exploring the role of environmental and anthropogenic changes[END_REF]. Shape analyses were performed using the free software R (R Development Core Team 2018). Fourier analysis is an efficient method used to describe teeth morphology [START_REF] Renaud | Size and shape variability in relation to species differences and climatic gradients in the African rodent Oenomys[END_REF].

The registration process proposed in this article relies on an algorithm of functional generalized Procrustes analysis (FGPA), which is an extension of the GPA algorithm that can be found in Dryden and Mardia (1998: p. 90) but is adapted in cases where the outline of an object is considered a continuous closed curve. The idea behind this modification is to construct a landmark-free registration method avoiding the delicate choice of the number and the position of landmarks. Originally, tooth outlines were described by a different number of points depending on the specimen (between 1500 and 2000). Then, the contour detection algorithm performed an optimized resampling of 250 points. The final number of points was arbitrarily chosen to be the most informative while limiting the dimensionality of the data. Consider an observed closed contour of a tooth arriving as M pairwise coordinates (x 1 , y 1 ),…, (x M , y M ) . We first considered that a contour line could be expressed as two curves (X(t), Y(t)) that are supposed to be a linear combination of K known basis functions such that:

where αk s and βk s are coefficients that must be estimated with regression on the sampled data, and K is the number of coefficients. The time index t arbitrarily belongs to [0;2π].

The basis functions ϕk, k = 1,⋯ , K are chosen by the practitioner and constitute a Fourier basis in our case. Other choices (B-splines, polynomials, …) can be relevant as well. The regression procedure was the same as that used in classical regression but in the functional case (see [START_REF] Ramsay | Functional Data Analysis[END_REF] for more details) moving from sampled contours to continuous curves. The number K of basis coefficients was arbitrarily fixed in such a way that a sufficient amount of curve variability was captured. In our case, 30 basis coefficients harmonics) were sufficient to describe the contour of the teeth, and additional coefficients did not provide more information.

Once every contour line was expressed through its estimated Fourier coefficients, it was possible to apply the registration FGPA algorithm using these coefficients as inputs. Registration of closed curves involves translating, rescaling, rotating and changing the phase of the configurations relative to each other to minimize an objective function (a total sum of squares). The registration of the starting points (phase changing) between curves was the main difference compared to the original Procrustes algorithm. They were not arbitrarily set and were part of the optimization process when registering the curves with Procrustes analysis, as shown in Fig. 3. See Online Resource 2 for more details about the registration algorithm.

Once the contour lines were registered, classical shape analyses could be achieved on the registered Fourier coefficients. To account for the entire shape variability, we used all Fourier coefficients in the following statistical analyses.

All subsequent treatments were applied to the two datasets (m1 and M1) separately, and all species were analysed together using the same protocol. To reduce data dimensionality in further analyses, we performed principal component analysis on both datasets using the PCA function of the FactoMineR package [START_REF] Le | FactoMineR: an R Package for multivariate analysis[END_REF]. We retained 95% of the shape variability, and the latest principal components are usually considered neglectable [START_REF] Baylac | Fourier descriptors, procrustes superimposition, and data dimensionality: an example of cranial shape analysis in modern human populations[END_REF]. In total, we used 20 principal components for the m1 dataset and 20 for the M1 dataset. As an estimate of tooth size, we used the surface area of masks.

Allometry

Allometry is the part of the shape caused by the influence of size [START_REF] Gould | Allometry and size in ontogeny and phylogeny[END_REF][START_REF] Klingenberg | Multivariate Allometry[END_REF][START_REF] Klingenberg | Size, shape, and form: concepts of allometry in geometric morphometrics[END_REF]. Size variation is an important determinant of morphological changes and consequently plays a critical role in the evolution of shape [START_REF] Klingenberg | Size, shape, and form: concepts of allometry in geometric morphometrics[END_REF]. If the FGPA procedure separates shape from size, it does not remove allometry. The presence of allometry was tested by performing the multivariate regression of shape on size using the function procD.lm from the library geomorph [START_REF] Adams | geomorph: an R package for the collection and analysis of geometric morphometric shape data[END_REF].

Clustering Analyses

To partition datasets without any prior information, we performed unsupervised clustering analyses based on shape variables. We chose to explore two complementary clustering protocols (summarized in Fig. 4):

1. Method 1 partitioned the global variability into MGs and gave an overview of MG evolution over time by retaining identity between MGs from one layer to another. However, this method implies that the phenotypic evolution that occurred in MGs over the time interval considered is less important than the phenotypic difference between MGs. Clustering was applied to the overall datasets.

2. Method 2 partitioned variability on each layer separately and thus considered phenotypic evolution occurring over the considered time interval. In contrast to method 1, method 2 is not based on any assumption about phenotypic evolution in MGs. However, it does not allow us to follow MG evolution over time because the relationship of MG identity determined in different layers is not known. The clustering was applied independently on each layer.

Clustering analyses were performed using the unsupervised morphological K-nearest neighbours method (KNN). KNN is a nonparametric classification (i.e., a method that is not based on statistical distributions) and therefore was applicable to both of our approaches (each of them implying highly varied sample sizes). The KNN classification is based on the assumption: "tell me who your neighbours are, and I will tell you who you are". In other words, each shape object is affected by its nearest neighbours' clusters. Unsupervised KNN was implemented in the clues function of the clues package [START_REF] Wang | CLUES: a non-parametric clustering method based on local shrinking[END_REF]). The dissimilarity measure used was Euclidean. In clues, the number of clusters and K (the number of neighbours to consider) were estimated by the algorithm. The number of clusters was obtained thanks to a partition procedure preceded by a local shrinking procedure, in which data points were "shrunk" towards a cluster centre [START_REF] Wang | CLUES: a non-parametric clustering method based on local shrinking[END_REF]). Then, K was selected between 1 and n-1 (n being the number of objects) based on the more robust clustering result [START_REF] Wang | CLUES: a non-parametric clustering method based on local shrinking[END_REF]. This robustness was assessed by the Silhouette index (SI) [START_REF] Kaufman | Finding Groups in Data: An Introduction to Cluster Analysis[END_REF], which measured the strengths of clusters. The SI was determined for each data point. It is comprised between -1 and 1. If SI > 0, the data point is closer to its assigned cluster than to other clusters. If SI < 0, the data point is misassigned to its cluster. If SI = 0, the data point is at equal distance from its assigned cluster and neighbouring clusters. To ensure that points are correctly allocated to clusters, the average SI was calculated and must be strictly positive [START_REF] Wang | CLUES: a non-parametric clustering method based on local shrinking[END_REF].

Clustering results and species are presented on PCA per layer in Online Resources 3, 4, 5 and 6. Size differences between MGs were tested using a one-factor ANOVA (analysis of variance) performed with the aov function of the stats library.

Shape Indicators

We studied three shape indicators over time: number of MGs, disparity and mean shape. 1. The number of MGs is the number of significantly distinct morphological units in a dataset. This value was obtained directly from the clustering analysis.

2. Disparity is the range of morphological variety in a group. It indicates changes in the amount of variation of a morphotype [START_REF] Gould | Wonderful life: the Burgess Shale and the Nature of History[END_REF][START_REF] Gould | The disparity of the Burgess Shale arthropod fauna and the limits of cladistic analysis: why we must strive to quantify morphospace[END_REF][START_REF] Wills | Disparity as an evolutionary index: a comparison of Cambrian and Recent arthropods[END_REF][START_REF] Foote | The evolution of morphological diversity[END_REF]). When considered independently of phylogenetic relationships, a good estimate of disparity is to quantify the amount of occupied morphospace [START_REF] Wills | Disparity as an evolutionary index: a comparison of Cambrian and Recent arthropods[END_REF]. As a disparity estimator, we used the median of the distances from the centroid of the group, a variance measure relatively insensitive to outliers [START_REF] Guillerme | Disparities in the analysis of morphological disparity[END_REF]. Disparity was obtained using the function dispRity.per.group from the package DispRity [START_REF] Guillerme | dispRity: a modular R package for measuring disparity[END_REF]. This function contains a bootstrap procedure to reduce the effect of sample composition, a common issue when studying disparity [START_REF] Butler | How do geological sampling biases affect studies of morphological evolution in deep time? A case study of pterosaur (Reptilia: Archosauria) disparity: geological sampling biases and disparity[END_REF]. Regarding sample size, variance is generally not biased by it [START_REF] Foote | The evolution of morphological diversity[END_REF], and no sample size correction was applied. We tested the statistical significance of disparity variation between the types of environments and between EH2 layers. Regarding method 1, we performed Friedman tests. This test is a rank sum test often used as an alternative to one-factor ANOVA on paired samples when the normality assumptions are not met. We used the friedman. test function from the stats library. Regarding method 2, we performed Kruskal-Wallis tests, a test similar to the Friedman test but adapted to independent samples. We used the Kruskal.test function from the stats library.

3. Mean shape reflects changes in the global shape of a group over time. It can detect shape divergence in groups in some layers. To quantify shape variations, we built a distance tree between mean shapes per MG per layer. We computed Mahalanobis distances and built the tree using the neighbour-joining algorithm, an agglomerative clustering method. The distance between groups in the tree was proportional to the morphological differences. To do that, we used the functions dist from the package stats and nj from the package ape [START_REF] Paradis | APE: Analyses of Phylogenetics and Evolution in R language[END_REF]). Visualizations of mean shapes per MG were also obtained. The neighbour joining distance trees were computed only for MGs obtained with method 1, as the MGs of method 2 had no continuity between layers. For each MG, we also tested the statistical significance of mean shape variation between the types of environments and between EH2 layers. We performed MANOVA (multivariate analysis of variance) and conducted pairwise Hotelling's T2 tests (the multivariate alternative to the t test). To correct the Hotelling's T2 results for multiple testing, we applied Bonferroni's criteria, which consists of identifying the number of tests n and then performing each of the n tests not at the α level of significance but at the α/n level. We used the manova function of the stats library and hotelling.test function from the Hotelling library (Curran 2017).

Biodiversity Indices

Biodiversity indices are frequently used in ecology, and some of them, such as the Shannon and Simpson indices, are also often applied to rodent fossil assemblages [START_REF] Avery | Micromammals as paleoenvironmental indicators and an interpretation of the late Quaternary in the Southern Cape Province[END_REF][START_REF] Geraads | The rodents from the late middle Pleistocene hominidbearing site of J'bel Irhoud, Morocco, and their chronological and paleoenvironmental implications[END_REF][START_REF] Meunier | Mise a jour de la liste des rongeurs d'Algerie, biogeographie et implications paleoecologiques[END_REF]. Those indices compute quantitative parameters about the structure of the specific diversity of a faunal assemblage (such as species richness, diversity and specific dominance). We choose three diversity indicators:

1. Number of taxa over time.

2. The Shannon-Weaver index (H'), which allows the estimation of the diversity of a community of organisms. A value close to 0 indicates that only one species is represented or that most organisms belong to the same species. The index is highest when many species are represented and when organisms are well distributed among species [START_REF] Blondel | Avifaune et vegetation. Essai d'analyse de la diversite[END_REF][START_REF] Hill | Diversity and evenness: a unifying notation and its consequences[END_REF]). This index is sensitive to variations in abundance of rare species [START_REF] Peet | The measurement of species diversity[END_REF] and highly sensitive to sample size when it is composed of fewer than 25 individuals (Cruz-Uribe 1988) (which was not the case here).

3. The Simpson index (D) measures the probability that two randomly selected individuals are of the same species. We used the unbiased Simpson index (D'), which is corrected for sample size, as suggested in [START_REF] Faith | The measurement of taxonomic evenness in zooarchaeology[END_REF]. It ranges from 0 (when diversity is maximized) to 1 (when diversity is minimized). The transformation 1-D' provides values that are easier to manipulate and interpret [START_REF] Pielou | An Introduction to Mathematical Ecology[END_REF][START_REF] Hill | Diversity and evenness: a unifying notation and its consequences[END_REF]). This index is sensitive to variations in abundance in more important species [START_REF] Peet | The measurement of species diversity[END_REF]. Biodiversity indices were computed for Rodentia alone and for all terrestrial microvertebrates, including Rodentia, Squamata, Amphibia, Chiroptera, Soricomorpha, and Erinaceomorpha, based on data from [START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF] and [START_REF] Stoetzel | Taphonomy and palaeoecology of the late Pleistocene to middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF]Stoetzel et al. ( , 2012b)).

Results

Allometry

A significant, albeit weak, allometric signal was found in both rodent m1 (R2 = 0.09; p value = 0.001**) and rodent M1 (R2 = 0.10; p value = 0.001**).

Clustering Analyses

The clustering results for both methods and both datasets are summarized in Table 2. All values of the SI means were in the range of (0,1), meaning data points were rightly assigned to their clusters. The detailed composition of each MG is available in Online Resources 7, 8, 9 and 10.

With method 1, we identified three MGs in the m1 dataset. MG1 contained mostly Mus spretus teeth but also some Dipodillus campestris, Meriones shawii/grandis and Apodemus sylvaticus teeth. MG2 and MG3 contained mostly Meriones shawii/grandis teeth but also some Dipodillus campestris and Mus spretus teeth. In the M1 dataset, two MGs were identified. MG1 was transgeneric and contained Mus spretus and Meriones shawii/grandis teeth, while MG2 was mainly composed of Meriones shawii/grandis teeth. MGs displayed significantly different sizes for m1 (ANOVA: p value < 2e-16 ***; df = 2; F = 355) and M1 (ANOVA: p value < 2e-16 ***; df = 1; F = 70.57).

For method 2, no identity between MGs from different layers was established. We identified up to three MGs per layer in the m1 dataset. Except for L7, all layers displayed only two MGs. Half of the MGs were transgeneric, while the other half were mainly composed of Meriones shawii/grandis teeth. In the M1 dataset, up to five MGs were identified. For m1, half of the MGs were transgeneric, while the other half were mainly composed of Meriones shawii/grandis teeth. MGs displayed significantly different sizes for m1 (ANOVA: p value < 2e-16 ***; df = 1; F = 145.9) and M1 (ANOVA: p value < 2e-16 ***; df = 1; F = 95.82).

MG size over time was explored and found to be invariant and thus uninformative with respect to our research question. Consequently, those results are neither presented nor discussed. For additional information, boxplots of lower and upper molar sizes per species are available in Online Resource 11.

Shape Indicators

The number of MGs and the disparity obtained for m1 over the eight layers of EH2 are presented in Fig. 5, and those obtained for M1 are presented in Fig. 6. Mean shape distance trees are presented in Fig. 7. For clarity, each analysis result is referred to as the "dataset-method", as in m1-method1.

The results showed several significant phenotypic variations. Global disparity variation over EH2 layers was statistically significant for m1-method1 (Friedman test: p value = 0.017*, df = 7, X2 = 17) but not for M1-method1 (Friedman test: p value = 0.19, df = 7, X2 = 10), m1-method2 (Kruskal-Wallis test: p value = 0.21, df = 7, X2 = 9.66) or M1-method2 (Kruskal-Wallis test: p value = 0.78, df = 7, X2 = 3.97). In m1-method1, we observed a global increase in disparity from L8 to L4a, followed by an important decline between L4a and L3 and then a slight increase until L1 was shared by all MGs (Fig. 5b). If we focused on MG disparity, this pattern of variation was also present in one/several MGs of m1-method2 (Fig. 5d) and in MG2 of M1-method1 (Fig. 6b). The variation in global disparity between environmental types was not statistically significant for m1-method1 (Friedman test: p value = 0.45, df = 3, X2 = 2.6), M1-method1 (Friedman test: p value = 0.14, df = 3, X2 = 5.4), m1-method2 (Kruskal-Wallis test: p value = 0.21, df = 3, X2 = 4.6) or M1-method2 (Kruskal-Wallis test: p value = 0.98, df = 3, X2 = 0.19).

Concerning the number of MGs over time, there was an increase in m1-method 2 and M1-method 2 in L7 (Figs. 5c, and6c), the only layer characterized by an open and arid environment. In M1method 2, there was also an additional MG in L5 (Fig. 6c), a layer characterized by an arid and semiopen environment.

Regarding mean shapes, the MGs of method 1 did not display the same internal shape variability. In m1-method1, MG1 showed more similar mean shapes, while MG3 had more diverse mean shapes (Fig. 7a). Similarly, in M1-method1, MG2 was composed of more similar mean shapes (except for the one in L1) than MG1 (Fig. 7b). In L3, the MGs of M1-method1 showed more similar mean shapes than those in other layers (Fig. 7b). Over the EH2 layers, the mean shape variation within each MG was statistically significant for m1-method1 (MANOVA: MG1: p value = 1.68e-9 ***, df = 7, F = 1.91, Pillai = 0.53; MG2: p value = 1.94e-8 ***, df = 7, F = 1.84, Pillai = 0.60; MG3: p value = 4.31e-06 ***, df = 7, F = 1.67, Pillai = 0.95) and for M1-method1 (MANOVA: MG1: p value = 4.76e-8 ***, df = 7, F = 1.81, Pillai = 0.68; MG2: p value < 2.2e-16 ***, df = 7, F = 3.069, Pillai = 0.64). The results of pairwise Hotelling's T2 tests are available in Online Resource 12. Variations between environmental types were also statistically significant for m1-method1 (MANOVA: MG1: p value = 1.61e-09 ***, df = 8, F = 2.56, Pillai = 0.30; MG2: p value = 0.00058 ***, df = 8, F = 1.73, Pillai = 0.24; MG3: p value = 0.00093 ***, df = 8, F = 1.72, Pillai = 0.42) and for M1-method1 (MANOVA: MG1: p value = 0.0023 **, df = 8, F = 1.63, Pillai = 0.26; MG2: p value < 2.2e-16 ***, df = 8, F = 4.30, Pillai = 0.37). The mean shape differences between environments per MG showed that the different MGs were not sensitive to the same environmental differences (Fig. 7). The results of pairwise Hotelling's T2 tests are presented in Online Resource 13.

Biodiversity Indices

Biodiversity indices over the eight studied stratigraphic layers of EH2 are presented in Fig. 8. All indicators displayed a continuous increase in diversity from L8 until L5-L4a, then a decrease until L2 and finally a second increase in L1 (Fig. 8). This general trend was nuanced in L7 by a decrease in the number of taxa of terrestrial vertebrates (Fig. 8A) and in L6 by a decrease in rodent diversity indicated by the Shannon and Simpson indices (Fig. 8b andc).

Discussion

We aimed to assess short-term multiscale phenotypic evolution in EH2 cave rodents relative to environmental changes. To do so, we studied the shape variation of rodents' teeth over time through an approach that set MGs instead of taxa as the phenotypic units. First, we will discuss the evidenced phenotypic variation over time. Then, we will discuss our methodological approach.

Phenotypic Variations Under Changing Environments

The main disparity variation consisted of an increase from L8 to L4a followed by a very large diversity drop between L4a and L3 and then an increase until L1 (around the late Pleistocene-Holocene transition). This result is illustrated by the disparity variations of the three MGs of m1-method1 and of MG2 of M1-method1, meaning there were important changes in the amount of variation within concerned morphotypes. The Shannon index of terrestrial microvertebrates displayed similar variation over time. This result supports the idea that the morphological variability in rodent teeth is consistent with variations in terrestrial microvertebrate diversity. Moreover, it demonstrated that the exclusion of some rare species from our sampling (Jaculus cf. orientalis and Eliomys cf. munbianus) did not significantly affect our results.

The ecomorphological meaning of those variations in disparity depends on the ecological and functional significance of the studied morphological features [START_REF] Ricklefs | Ecological and evolutionary inferences Reilly SM[END_REF][START_REF] Van Valkenburgh | Ecomorphological analysis of fossil vertebrates and their paleocommunities[END_REF]. In the case of teeth, because morphological changes can be related to the environment through diet [START_REF] Wolf | Does weather shape rodents? Climate related changes in morphology of two heteromyid species[END_REF][START_REF] Mcguire | Geometric morphometrics of vole (Microtus californicus) dentition as a new paleoclimate proxy: shape change along geographic and climatic clines[END_REF][START_REF] Coillot | Correlated changes in occlusal pattern and diet in stem murinae during the onset of the radiation of old world rats and mice[END_REF][START_REF] Cano | Dietary ecology of murinae (Muridae, Rodentia): a geometric morphometric approach[END_REF][START_REF] Rodrigues | The great variety of dental structures and dynamics in rodents: new insights into their ecological diversity[END_REF][START_REF] Pineda-Munoz | Inferring diet from dental morphology in terrestrial mammals[END_REF], changes in disparity/diversity patterns can reflect environmental variations [START_REF] Erwin | Disparity: morphological pattern and developmental context[END_REF]. The fact that disparity variations in our results were displayed by two morphological features (m1 and M1) and were shared by MGs that included diverse species (Meriones shawii/grandis, Dipodillus campestris and Mus spretus) suggested that it was driven by (an) external disturbance(s). However, it is difficult to associate specific phenotypic responses to particular environments. Moreover, no major environmental change was identified by paleoenvironmental studies between L4a and L3, while an important variation in disparity was identified. Moreover, clades other than rodents were impacted: variations in shrew species composition have been observed in the same layers [START_REF] Cornette | Shrews of the genus Crocidura from El Harhoura 2 (Temara, Morocco): the contribution of broken specimens to the understanding of Late Pleistocene-Holocene palaeoenvironments in North Africa[END_REF]. Thus, the considered environmental conditions do not seem to be entirely related to disparity variations. Either disturbances of another nature-ecological or anthropogenic-must have caused those phenotypic changes, or there have been environmental changes that could not be detected through the usual markers used in paleoenvironmental inferences. The disparity variation of m1-method1 may support this last hypothesis, as it was statistically significant between layers but not between environmental types.

This mismatch might be explained by the fact that paleoenvironmental inferences can sometimes be questionable. At EH2, paleoclimatic fluctuation assessments were based on standard analyses of microvertebrate assemblages and were complemented by stable carbon (13C) and oxygen (18O) isotope analyses on Meriones shawii/grandis teeth [START_REF] Jeffrey | Exploring palaeoaridity using stable oxygen and carbon isotopes in small mammal teeth: a case study from two Late Pleistocene archaeological cave sites in Morocco, North Africa[END_REF]. These analyses indicated that L11 to L7 (dated to MIS 5) encountered more humid conditions than today, while L6 to L3 (dating from MIS 5 to 3) were subjected to conditions similar to the present conditions. This result suggests a relatively open environment throughout the Late Pleistocene. L2 (dated to MIS 2 and/or beginning of MIS 1) appeared more arid, which was consistent with the aridification of the region at this period. We thus observed some discrepancies with the previous palaeoecological data, notably in L7, L6 and L5. Isotopic results indicated relatively humid conditions in L7 and L5 and more arid conditions in L6, while faunal communities indicated the opposite. These discrepancies could be explained by different phenomena. One of these could be the differences in the signals recorded from a set of species communities and those recorded from individuals within a single species (which can adapt differently to environmentalchanges) [START_REF] Stoetzel | Quaternary terrestrial microvertebrates from Mediterranean northwestern Africa: state-of-the-art focused on recent multidisciplinary studies[END_REF]. Consequently, the paleoenvironmental signal at the individual or species/MG scales may be different than that at the community scale.

An alternative explanation for this mismatch might reside in the fact that dental morphology reflects the environment through diet [START_REF] Wolf | Does weather shape rodents? Climate related changes in morphology of two heteromyid species[END_REF][START_REF] Mcguire | Geometric morphometrics of vole (Microtus californicus) dentition as a new paleoclimate proxy: shape change along geographic and climatic clines[END_REF][START_REF] Coillot | Correlated changes in occlusal pattern and diet in stem murinae during the onset of the radiation of old world rats and mice[END_REF][START_REF] Cano | Dietary ecology of murinae (Muridae, Rodentia): a geometric morphometric approach[END_REF][START_REF] Rodrigues | The great variety of dental structures and dynamics in rodents: new insights into their ecological diversity[END_REF][START_REF] Selig | Mammalian molar complexity follows simple, predictable patterns[END_REF]. However, the dietary preferences of the species considered may not be directly linked to environmental inferences. Some rodent species in our sample have flexible diets (e.g., [START_REF] Palomo | Mus spretus (Rodentia: Muridae)[END_REF][START_REF] Adamou-Djerbaoui | Etude du regime alimentaire d'un rongeur nuisible (Meriones shawii Duvernoy, 1842, Mammalia, Rodentia) en Algerie[END_REF] and thus can display broad dietary categories. Then, their diet composition depends not only on resources available but also on other ecological factors, such as competition [START_REF] Ledevin | Phylogeny and adaptation shape the teeth of insular mice[END_REF] indicates that competition with other species may be a more important driver of insular mouse molar shape than the environment. Consequently, the shape indicators computed on MGs and environmental inferences may not be related to the exact same environmental and ecological parameters. Moreover, MGs may not be exclusively partitioned on ecological characteristics; rather, they may also be partitioned on life history traits, meaning they may not display equal phylogenetic signals. For example, MG2 of M1-method1 was mainly composed of one species complex, Meriones shawii/grandis, suggesting that it might have been clustered based on traits that reflect history. This would explain why MG was not correlated with environmental inferences.

Regarding the m1 dataset, the MG comprising Mus spretus (MG1) displayed significantly similar mean shapes per layer (except in L1). This result suggests a high morphological stability of this MG over time and with environmental changes. This result is consistent with the continuous presence of this species over the EH2 sequence, with no significant changes in terms of abundance over time [START_REF] Stoetzel | Taphonomy and palaeoecology of the late Pleistocene to middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF]) and a stability in morphology and size of molars from the Late Pleistocene to the Middle Holocene [START_REF] Stoetzel | Mus in Morocco: a Quaternary sequence of intraspecific evolution[END_REF]. This result may be related to a strong adaptability of this generalist species to environmental changes and/or a continuous presence of suitable habitats for this species in the Rabat-Temara region through the Late Quaternary [START_REF] Lalis | Out of Africa: demographic and colonization history of the Algerian mouse (Mus spretus Lataste)[END_REF]. Geometric morphometric analysis of Meriones molars has shown that during the Late Pleistocene, both M. shawii and M. grandis were present in western Morocco [START_REF] Stoetzel | Systematics and evolution of the Meriones shawii/grandis complex (Rodentia, Gerbillinae) during the Late Quaternary in northwestern Africa: exploring the role of environmental and anthropogenic changes[END_REF]. Few changes occurred throughout this period despite several climatic changes, indicating that suitable habitats probably always occurred in the region, as was the case for Mus spretus. However, during the Middle Holocene, paleontological and morphometric analyses highlighted a change in the faunal spectrum and species distribution accompanied by a clear demographic collapse [START_REF] Lalis | Phylogeography and demographic history of Shaw's Jird (Meriones shawii complex) in North Africa[END_REF][START_REF] Stoetzel | Adaptations and dispersals of anatomically modern humans in the changing environments of North Africa: the contribution of microvertebrates[END_REF]. The authors therefore suggest that during the humid period of the Middle Holocene, when Mediterranean forests expanded, the semiarid open steppe habitats, which were more suitable for Meriones, decreased drastically in the northern half of Morocco, leading to population collapses and changes in their geographic distribution. Consequently, we hypothesize that Mus spretus may have been less affected by environmental changes than Meriones shawii/grandis and that the identified disturbance(s) in L4a must be of a climatic nature. This result is supported by isotope studies that display higher δ13C values than the normal values observed at EH2, possibly indicating an important change in vegetation at this time [START_REF] Jeffrey | Exploring palaeoaridity using stable oxygen and carbon isotopes in small mammal teeth: a case study from two Late Pleistocene archaeological cave sites in Morocco, North Africa[END_REF][START_REF] Stoetzel | Quaternary terrestrial microvertebrates from Mediterranean northwestern Africa: state-of-the-art focused on recent multidisciplinary studies[END_REF]. Another possibility would be to relate this change in fauna with predator activity. The presence of variable nonhuman predators (especially small mammalian carnivores and birds of prey) has been attested at EH2 cave and is suspected to be the main agent of modification for faunal remains [START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF][START_REF] Campmas | Were Upper Pleistocene human/nonhuman predator occupations at the Temara caves (El Harhoura 2 and El Mnasra, Morocco) influenced by climate change?[END_REF][START_REF] Campmas | Which predators are responsible for faunal accumulations at the Late Pleistocene layers of El Harhoura 2 Cave (Temara, Morocco)?[END_REF]. Even if the microvertebrate diversity does not appear to be related to predator shifts between levels [START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF][START_REF] Stoetzel | Taphonomy and palaeoecology of the late Pleistocene to middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF], microevolutionary trends might have been affected, thereby influencing macroevolutionary trends [START_REF] Hansen | Translating between microevolutionary process and macroevolutionary patterns: the correlation structure of interspecific data[END_REF].

The combination of the different shape indicators can allow us to formulate hypotheses about evolution scenarios. This was the case for the phenotypic changes in M1-method1 that occurred in L3, and especially MG2. In this layer, there is a consistent decline in disparity, combined with a mean shape convergence between the two MGs. Mean shapes represent clumps of forms in the occupied morphospace and disparity the dispersal of those clumps. It is known that morphological clumps can represent local optima surrounded by less fit alternatives [START_REF] Erwin | Disparity: morphological pattern and developmental context[END_REF]. Shifts in mean shapes should then happen when there is a shift in ecology and thus in optima. Moreover, changes in disparity patterns can reveal patterns of morphological selectivity [START_REF] Erwin | Disparity: morphological pattern and developmental context[END_REF]. Thus, a reduction in disparity combined with a mean shape shift towards a particular morphotype can be interpreted as the selection of a form that is advantageous with respect to environmental conditions. Then, we can hypothesize that a morphological selective event occurred in L3, resulting in a shape convergence between the two MGs.

In L7 in m1-method 2 and M1-method 2, an increase in the number of MGs coupled with a decrease in the disparity of some MGs was observed. The presence of a new morphotype might have reinforced competition for resources, resulting in each morphotype converging towards its mean shape, reducing each group's disparity. Faunal variations in this layer were also observed in shrews with a particular species composition in the layer [START_REF] Cornette | Shrews of the genus Crocidura from El Harhoura 2 (Temara, Morocco): the contribution of broken specimens to the understanding of Late Pleistocene-Holocene palaeoenvironments in North Africa[END_REF]. This result might be explained by a particularly arid and open environment at this time [START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF][START_REF] El Hajraoui | Prehistoire de La Region de Rabat-Temara[END_REF]. However, isotope studies disagree with this hypothesis and suggest there were more humid conditions [START_REF] Jeffrey | Exploring palaeoaridity using stable oxygen and carbon isotopes in small mammal teeth: a case study from two Late Pleistocene archaeological cave sites in Morocco, North Africa[END_REF]. The hypothesis of an arid climate in L7, as in L5, is mainly supported by the species abundance and the presence of Jaculus cf. orientalis. However, this steppic species can also be considered an indicator of more continental conditions (i.e., shoreline retreat caused by a sealevel drop) rather than of particularly arid conditions. Consequently, some arguments question the paleoenvironmental inferences in L7 based on faunal communities, even if it is almost undeniable that an environmental change occurred.

In L6, MG1 of M1-method1 displayed a very large decline in disparity. An important fact to underline is that L7 and L6 (as for L4a and L3) are very close layers in the data [START_REF] Jacobs | Single-grain OSL chronologies for Middle Palaeolithic deposits at El Mnasra and El Harhoura 2, Morocco: implications for Late Pleistocene human-environment interactions along the Atlantic Coast of northwest Africa[END_REF][START_REF] Janati-Idrissi | Datation par ESR-U/th combinees de dents fossiles des grottes d'El Mnasra et d'El Harhoura 2, region de Rabat-Temara. Implications chronologiques sur le peuplement du Maroc atlantique au Pleistocene superieur et son[END_REF]Ben Arous et al. 2020a, 2020b). Thus, these phenotypic variations have a unique cause. We must also keep in mind that two species, which we know are present in these deep layers, were not included in this study (Jaculus cf. orientalis in L5 and L7 and Eliomys cf. munbyanus in L5 and L8). Therefore, the importance of these trends might be underestimated by this bias in the sampling.

Lower and Upper Molars, Two Complementary Models

The shapes of the two biological objects-m1 and M1-globally reflected evolutionary patterns similar to those observed in other rodent species (e.g., [START_REF] Renaud | Morphological evolution of the murine rodent Paraethomys in response to climatic variations (Mio-Pleistocene of North Africa)[END_REF]. More diversity was identified among m1, but more varied trends were displayed by M1. Statistical tests of the mean shape differences between environmental types indeed showed that M1 was more sensitive to environmental changes than was m1. However, these results were eventually less informative because they were less specific to particular environmental transitions. Thus, we did not find a clear higher responsiveness of M1 to environmental changes, contrary to what was concluded by [START_REF] Renaud | Differential evolvability along lines of least resistance of upper and lower molars in island house mice[END_REF]. However, this result makes sense, as [START_REF] Renaud | Differential evolvability along lines of least resistance of upper and lower molars in island house mice[END_REF] explained this stronger responsiveness was due to the fact that M1 has a higher allometric component than m1, inducing a greater evolvability of M1, and that in our dataset, both M1 and m1 display similarly weak allometry. Indeed, we found that tooth size explained only ~ 10% of shape variation in both m1 and M1, which was quite low but was within the range of allometry usually found in rodent teeth (e.g., ~ 3% in insular mice, as found by [START_REF] Ledevin | Phylogeny and adaptation shape the teeth of insular mice[END_REF], and ~ 20% in the Orkney vole, as found by [START_REF] Cucchi | The changing pace of insular life: 5000 years of microevolution in the orkney vole (Microtus arvalis orcadensis))[END_REF]). Each of the biological models allowed us to detect different particular trends: different phenotypic responses within one species in the case of m1 and, for M1, a morphological selective event in L3. These two phenotypic proxies thus provide complementary information leading to a better understanding of rodent phenotypic evolution.

Biological Interpretation of the MGs

The MGs are taxon-free groups of specimens with more similar teeth shapes and significantly different sizes. They were also characterized by phenotypic evolutionary characteristics. For example, the MGs of m1-method1 displayed varied mean shape similarities between layers. This could be interpreted as different degrees of morphological stability over time and through varied environments. MG1 had significantly similar mean shapes (except in L1), possibly implying a high resistance to external variations. Conversely, MG2 had more significantly different mean shapes between layers, meaning that this MG was more plastic. Moreover, MGs were not equally sensitive to all environmental changes, suggesting that they were responsive to different environmental parameters. MG1 of m1-method1 was sensitive only to transitions between humid and close environments to all other types of environments. MG1 of M1-method1 was very insensitive to changes in the environment, while MG2 was highly sensitive. Thus, MGs seemed to have different responses to disturbances. They might be interpreted as taxon-free ecological response groups, which can be expected because rodents are ecological opportunists with flexible diets and habitat preferences [START_REF] Nowak | Walker's Mammals of the World[END_REF][START_REF] Alhajeri | Ecological and ecomorphological specialization are not associated with diversification rates in muroid rodents (Rodentia: Muroidea)[END_REF].

We detected more phenotypic variations than what is evidenced by conventional indicators of biodiversity, although our results remain consistent with the latter. This may be explained by the fact that during the Late Quaternary, environmental changes were less abrupt in North Africa than in Europe, and no major faunal turnover was observed. The geographical location of the EH2 cave also implies that fauna and vegetation are strongly influenced by the Atlantic ocean, which was probably the case during the whole Quaternary. This includes fewer variations in humidity and temperature on the coast than in the interior of the continent, as well as regional faunal differences [START_REF] Stoetzel | Adaptations and dispersals of anatomically modern humans in the changing environments of North Africa: the contribution of microvertebrates[END_REF][START_REF] Stoetzel | Quaternary terrestrial microvertebrates from Mediterranean northwestern Africa: state-of-the-art focused on recent multidisciplinary studies[END_REF]. In sum, only minor changes in microvertebrate communities were observed throughout the EH2 sequence [START_REF] Stoetzel | Taphonomy and palaeoecology of the late Pleistocene to middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Temara, Morocco)[END_REF][START_REF] Stoetzel | Late Cenozoic micromammal biochronology of northwestern Africa[END_REF][START_REF] Stoetzel | Adaptations and dispersals of anatomically modern humans in the changing environments of North Africa: the contribution of microvertebrates[END_REF]. This means that a 'common pool' of species was preserved through time with only a few differences in the composition and proportions of species. Microvertebrates experienced low-amplitude environmental changes, although alternations between arid and humid periods were recorded at EH2. The landscapes would have been characterized by a mosaic of habitats, with only the relative cover of the different habitats changing over time. Most of these changes have certainly been too subtle to be detected by the classic biodiversity indices.

In an archaeological context, this approach may represent another way of apprehending phenotypic diversity. It is recognized that in archaeology, global diversity is often underestimated because of the difficulty in identifying some taxa [START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF]. Moreover, even current taxa taxonomy can be unstable, as for the Meriones shawii/grandis complex, the systematic nature of which has been very controversial for a long time [START_REF] Carleton | Order Rodentia[END_REF]Darvish 2011;[START_REF] Lalis | Phylogeography and demographic history of Shaw's Jird (Meriones shawii complex) in North Africa[END_REF][START_REF] Stoetzel | Systematics and evolution of the Meriones shawii/grandis complex (Rodentia, Gerbillinae) during the Late Quaternary in northwestern Africa: exploring the role of environmental and anthropogenic changes[END_REF]. The MG approach is independent of identification, and in some cases, unsupervised clustering proved to better evaluate phenotypic diversity than species [START_REF] Quenu | Geometric morphometrics and machine learning challenge currently accepted species limits of the land snail Placostylus (Pulmonata: Bothriembryontidae) on the Isle of Pines, New Caledonia[END_REF]. Moreover, the association of geometric morphometrics and machine learning has repeatedly proven its ability in terms of phenotype discrimination [START_REF] Dubey | Potential of artificial neural networks in varietal identification using morphometry of wheat grains[END_REF][START_REF] Bocxlaer | Comparison of morphometric techniques for shapes with few homologous landmarks based on machine-learning approaches to biological discrimination[END_REF][START_REF] Cornette | Shrews of the genus Crocidura from El Harhoura 2 (Temara, Morocco): the contribution of broken specimens to the understanding of Late Pleistocene-Holocene palaeoenvironments in North Africa[END_REF][START_REF] Guillaud | Is vertebral form a valid species-specific indicator for salmonids? The discrimination rate of trout and Atlantic salmon from archaeological to modern times[END_REF][START_REF] Mapp | Otolith shape and size: the importance of age when determining indices for fish-stock separation[END_REF][START_REF] Soda | Artificial neural networks and geometric morphometric methods as a means for classification: a case-study using teeth from Carcharhinus sp. (Carcharhinidae)[END_REF][START_REF] Fang | Beak identification of four dominant octopus species in the East China Sea based on traditional measurements and geometric morphometrics[END_REF][START_REF] Quenu | Geometric morphometrics and machine learning challenge currently accepted species limits of the land snail Placostylus (Pulmonata: Bothriembryontidae) on the Isle of Pines, New Caledonia[END_REF]). As a consequence, it might be a complementary local diversity indicator to species, as it does not carry the same information and potentially represents ecological groups.

However, this approach has the disadvantage of being a relative rather than an absolute characterization of phenotypic diversity, which could sometimes lead to an underestimation of the overall diversity when compared with other localities [START_REF] Read | Biocrust morphogroups provide an effective and rapid assessment tool for drylands[END_REF]. Moreover, we must keep in mind that Jaculus cf. orientalis (present in L5 and L7) and Eliomys cf. munbyanus (present in L5 and L8) were not considered in the present study because of the lack of usable teeth for geometric morphometric analyses, which could have led to an underestimation of the diversity in the deepest levels.

The MG approach proposed in this paper has a second disadvantage: the sine qua non assumption of method 1. Herein, the MG approach was addressed through two complementary methods. The major difference between the two is that method 1 partitions the variability once on the overall dataset, while method 2 does it layer by layer. By pooling remains from different layers, the partitioning protocol of method 1 implies that phenotypic evolution occurring in MGs over time is less important than the phenotypic difference between MGs. It is questionable whether this hypothesis is fully satisfied. First, the number of MGs varied greatly between the two methods. Phenotypic evolution had already been brought to light over the studied time interval (i.e., between layers). Indeed, morphological differences were observed in remains of the same taxon between old and recent layers. These differences were attributed to the presence of "primitive" morphological characters in old layers [START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF]. However, the overall similarity of trends between the two methods suggests that it is doubtful that those phenotypic changes are equally or more important than phenotypic differences between MGs. The case of EH2 rodents seems to meet the hypothesis of method 1 to some extent. However, this was not the case for any dataset, and the MG approach will be applicable only to a limited number of cases.

Conclusion

We observed significant short-term shape changes within MGs at both intra-and interspecific scales and identified significant variations in MG disparities over time. Thus, MGs appeared to be phenotypic response units that could be transversal between species. However, environmental changes did not seem to be the main determinants of characterized phenotypic variations, as they were not related to changes in disparity. Nevertheless, discrepancies between paleoenvironmental proxies might suggest the existence of some uncertainties in paleoenvironmental inferences that might explain our results. Another lead may be to look at the functional characteristics of MGs rather than the morphological ones. Because they more directly represent the ability of individuals to perform ecologically relevant tasks, they might allow one to obtain a better characterization of the phenotypic responses to environmental changes. Supplementary Information. The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s10914-022-09620-w.

FIGURES

Fig. 1. Stratigraphy, age, and characteristic environments of the eight layers of EH2 used in this study, color-coded by environment [START_REF] Michel | La macrofaune du Pleistocene superieur d'El Harhoura 2 (Temara, Maroc) : donnees preliminaires[END_REF][START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF][START_REF] Jacobs | Single-grain OSL chronologies for Middle Palaeolithic deposits at El Mnasra and El Harhoura 2, Morocco: implications for Late Pleistocene human-environment interactions along the Atlantic Coast of northwest Africa[END_REF][START_REF] Janati-Idrissi | Datation par ESR-U/th combinees de dents fossiles des grottes d'El Mnasra et d'El Harhoura 2, region de Rabat-Temara. Implications chronologiques sur le peuplement du Maroc atlantique au Pleistocene superieur et son[END_REF][START_REF] Nespoulet | Excavation report[END_REF]Ben Arous et al. 2020a, b;[START_REF] Marquer | The first use of olives in Africa around 100,000 years ago[END_REF]. Three dating methods were used: AMS-14C based on organic remains, combined US-ESR (Combined uranium series and electron spin resonance) and OSL (optical stimulated luminescence). Layers unused in this study are in grey. Table 2. Synthesis of clustering results of method 1 (one general clustering performed) and method 2 (one clustering performed per layer). M1: upper first molars; m1: lower first molars.

Fig. 2 .

 2 Fig. 2. Acquisition protocol from picture taking to mask extraction. The protocol varied depending on whether the tooth was individualized or not (i.e., separate or still set into the mandible or the maxilla).

Fig. 3 .

 3 Fig. 3. Aligned outlines of right lower first molars with the registered starting point of each curve.

Fig. 4 .

 4 Fig. 4. Protocols of the two clustering methods.

Fig. 5 .

 5 Fig. 5. Number of MGs and disparity of m1 (lower first molar) shape over the eight studied stratigraphic layers of EH2. Environmental conditions (Stoetzel 2009; El Hajraoui et al. 2012) are indicated by background colors. a Number of MGs obtained with method 1; b disparity of MGs of method 1; c Number of MGs obtained with method 2; d disparity of MGs of method 2.

Fig. 6 .

 6 Fig. 6. Number of MGs and disparity of M1 (upper first molar) shape over the eight studied stratigraphic layers of EH2. Environmental conditions (Stoetzel 2009; El Hajraoui et al. 2012) are indicated by background colors. a Number of MGs obtained with method 1; b disparity of MGs of method 1; c Number of MGs obtained with method 2; d disparity of MGs of method 2.

Fig. 7 .

 7 Fig.7. Neighbor joining distance tree of the mean shape of MGs of method 1. Each layer of EH2 is referred as follows: "layer 1" is "L1". Environmental conditions[START_REF] Stoetzel | Les microvertebres du site d'occupation humaine d'El Harhoura 2 (Pleistocene superieur -Holocene, Maroc) : systematique, evolution, taphonomie et paleoecologie. Dissertation, Museum national[END_REF][START_REF] El Hajraoui | Prehistoire de La Region de Rabat-Temara[END_REF]) are indicated by background colors. a Neighbor joining distance tree between mean shape of MGs per layer for m1 (lower first molars) and results of pairwise Hotelling's T2 tests between environmental types per MGs; b Neighbor joining distance tree between mean shape of MGs per layer for M1 (upper first molars) and results of pairwise Hotelling's T2 tests between environmental types per MGs.

Fig. 8 .

 8 Fig. 8. Biodiversity indices over the eight studied stratigraphic layers of EH2. Environmental conditions (Stoetzel 2009; El Hajraoui et al. 2012) are indicated by background colors. a Number of taxa; b Shannon index; c Simpson index.

Acknowledgements. The authors would like to thank the Mission Archeologique El Harhoura-Temara, Ministere des Affaires Etrangeres et Europeennes, France, Ministere de la Culture, Morocco (dir. R. Nespoulet and M.A. El Hajraoui) and the Institut National des Sciences de l'Archeologie et du Patrimoine (INSAP), Rabat, Morocco for helping to provide access to EL Harhoura 2 fossil specimens. We thank Eslem Ben Arous for discussions about dating and chronological context of the El Harhoura 2 and neighboring sites. We also thank the two anonymous reviewers and Ana Rosa Gomez Cano for constructive comments that greatly helped to improve the manuscript.

Funding. This research was funded by the Universite 463 Paris Descartes, the Ecole Doctorale FIRE -Programme Bettencourt (L.T.). L.T. and E.S. received financial support from the ATM ArviClim project (E. Stoetzel & R. Cornette dir.) from the Museum national d'Histoire naturelle (MNHN), Paris, France, for a mission at the INSAP in Rabat in July 2019, allowing the study of additional specimens.

Availability of Data and Code. The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request. The R script used for analyses is available in Online Resource 14.

Declarations

Conflict of Interest The authors have no conflicts of interest to declare that are relevant to the content of this article.