Louis Arbez
email: alouis.arbez@u-bourgogne.fr.

Aurélien Royer

Danielle Schreve

Rémi Laffont

Serge David

Sophie Montuire

The missing Myopus: plugging the gaps in Late Pleistocene small mammal identification in western Europe with geometric morphometrics

Keywords: Geometric morphometrics, Lemming, Late Glacial, Paleobiogeography, Molar shape, Boreal environment, Taxonomy

Lemmus and Myopus are two lemming species with distinct habitat requirements but which show very similar dental morphologies. They are thus extremely difficult to distinguish from one another in the fossil record on the basis of their dental remains, leading to poor understanding of the palaeobiogeographical evolution of Myopus as well as inaccurate paleoenvrionmental reconstructions. Currently, the presence of Myopus in the fossil register from the Pleistocene is still debated and no firm occurrence of this lemming in western Europe has yet been confirmed during the Late Pleistocene. In this paper, we used geometric morphometrics on modern material to establish morphological differences between Lemmus and Myopus teeth (first lower and third upper molars). Morphological data was then used to build a robust linear discriminant model able to confidently classify isolated teeth of these two genera, and finally, linear discriminant models were used on fossil remains of Lemmus/Myopus from two Late Pleistocene archaeological/palaeontological sites (Grotte des Gorges and Gully Cave). This study demonstrates, for the first time, the presence of Myopus schisticolor in west European Late Pleistocene sites between the end of MIS 3 and the beginning of Holocene, during climatic events that favoured the development of taiga forest of birch and pine in these regions.

Introduction

Rodents are characterized by a significant taxonomic diversity associated with various ecological and biological traits, thereby ensuring their widespread occurrence throughout all modern biotopes. They have been intensively studied with respect to the fields of evolution, biochronology and palaeoenvironmental change (e.g. [START_REF] Chaline | Les rongeurs et l'évolution[END_REF][START_REF] Cox | Evolution of the Rodents: Advances in Phylogeny, Functional Morphology and Development[END_REF]Fejfar and Heinrich, 1990;Hernandez Fernandez, 2006;van Kolfschoten, 1995;[START_REF] Kowalski | Lemmings [Mammalia, Rodentia] as indicators of temperature and humidity in the European Quaternary[END_REF][START_REF] Montuire | Rodents and climate. 1. A model for estimating past temperatures using arvicolids (Mammalia: Rodentia)[END_REF][START_REF] Rekovets | Phenomenon in the Evolution of Voles (Mammalia, Rodentia, Arvicolidae)[END_REF]. Their remains, well preserved in many archaeological and palaeontological contexts, have therefore contributed to the reconstruction of patterns of Quaternary palaeoclimatic change, still a major challenge in our understanding of continental chronostratigraphies today. Among the methods applied for paleoclimate reconstructions, many are based on species associations, such as the presence/absence or relative abundance of species, hence the need to have accurate taxonomic identifications (e.g. [START_REF] Gobalet | A Critique of Faunal Analysis; Inconsistency among Experts in Blind Tests[END_REF][START_REF] Lyman | Taxonomic identification of zooarchaeological remains[END_REF][START_REF] Lyman | Assumptions and Protocol of the Taxonomic Identification of Faunal Remains in Zooarchaeology: a North American Perspective[END_REF][START_REF] Royer | New bioclimatic models for the quaternary palaearctic based on insectivore and rodent communities[END_REF][START_REF] Stahl | The recovery and interpretation of microvertebrate bone assemblages from archaeological contexts[END_REF].

Numerous obstacles hinder the identification of small-bodied animals, even when their skeletal remains are recovered [START_REF] Stahl | The recovery and interpretation of microvertebrate bone assemblages from archaeological contexts[END_REF]. One problem is that while modern zoologists have many criteria by which to identify animals (e.g. fur, size, skull and post-cranial anatomy, geographical range), most of these criteria cannot be applied to fossils. A second problem is that identification of small mammals is based on limited types of remains, generally the teeth, which are most resistant to mechanical taphonomical processes. Despite dental remains being the most diagnostic elements at the genus or species levels, teeth may exhibit very similar morphologies, which could potentially prevent confident identifications. The last problem also reflects the differential state of diagnostic information for relevant taxa [START_REF] Lyman | Assumptions and Protocol of the Taxonomic Identification of Faunal Remains in Zooarchaeology: a North American Perspective[END_REF][START_REF] Stahl | The recovery and interpretation of microvertebrate bone assemblages from archaeological contexts[END_REF]. Many criteria have been proposed to identify these smallest-size species, but most of them rely on the "holotypic" morphology (i.e. the typical form as defined by characteristic specimens of this species), without taking into account geographical or inter-population variability, except to point out very unusual cases. The ensuing large number of criteria proposed for species identification often leads either to disagreement or to differential use of these criteria among taxonomists. As said by [START_REF] Stahl | The recovery and interpretation of microvertebrate bone assemblages from archaeological contexts[END_REF], it must always be remembered that higher level taxonomic categories are necessarily dynamic and constantly open to revision.

Among the rodents, voles and lemmings, which show substantial tooth variability throughout species, time and space, are key taxa for tracking Quaternary environmental and climatic changes (e.g. [START_REF] Kowalski | Pleistocene rodents of Europe[END_REF][START_REF] Montuire | Molar shape differentiation during range expansions of the collared lemming (Dicrostonyx torquatus) related to past climate changes[END_REF][START_REF] Nadachowski | Late quaternary rodents of Poland with special reference to morphotype dentition analysis of voles[END_REF][START_REF] Royer | Etude paléoenvironnementale et paléoclimatique du Pléistocène supérieur du Sud-Ouest de la France, à partir d'analyses comparées d'associations fauniques et de biogéochimies effectuées sur les micromammifères[END_REF].

Lemmings are reliable biomarkers of specific and typical high latitude environments [START_REF] Chaline | Les rongeurs du Pléistocène moyen et supérieur de France[END_REF][START_REF] Kowalski | Lemmings [Mammalia, Rodentia] as indicators of temperature and humidity in the European Quaternary[END_REF]. In Eurasia, four lemming genera are found in fossil contexts: 1) Dicrostonyx and Lagurus, which are regularly present during the coldest or driest phases of the Late Pleistocene extending as far as the south of France [START_REF] Chaline | Les rongeurs du Pléistocène moyen et supérieur de France[END_REF]Marquet, 1993;[START_REF] Royer | Investigating the Influence of Climate Changes on Rodent Communities at a Regional-Scale (MIS 1-3, Southwestern France)[END_REF]; 2) Lemmus, which is also frequently found throughout Eurasia in open and cold environments, but still in modest abundance [START_REF] Kowalski | Lemmings [Mammalia, Rodentia] as indicators of temperature and humidity in the European Quaternary[END_REF][START_REF] Kowalski | Pleistocene rodents of Europe[END_REF], and 3) Myopus, which has so far only been attested in Russia at this time, associated with more temperate species (e.g. [START_REF] Markova | Institute of geography of Russian Academy of Sciences[END_REF]Ponomarev et al., 2013aPonomarev et al., , 2013bPonomarev et al., , 2013c)).

As Lemmus (Link, 1795) and Myopus (Miller, 1910) have two strong but divergent ecological signals, their taxonomic discrimination in the fossil assemblages is essential in order to investigate past faunal associations and understand their significance in paleoclimatic reconstructions. The True Lemming genus (Lemmus) is a complex of several species which have a northern distribution restricted to Arctic tundra environments associated with high seasonality and little vegetation. These taxa use the thick snow covering as a shelter, giving them access to the vegetation persisting under this natural shield during winter [START_REF] Domine | Snow physical properties may be a significant determinant of lemming population dynamics in the high Arctic[END_REF]Reid et al., 2011;[START_REF] Stenseth | Biology of lemmings[END_REF], as well as using bogs and wetlands during summer and eschewing forest as much as possible [START_REF] Vaillant | Spatial distribution in Norwegian lemming Lemmus lemmus in relation to the phase of the cycle[END_REF]Wilson et al., 2017).

The wood lemming (Myopus schisticolor [Lilljeborg, 1844]) is the only representative of its genus and inhabits exclusively coniferous forest with a thick moss cover, which constitutes its main food source [START_REF] Bobretsov | Population dynamics of wood lemming (Myopus schisticolor) in different landscapes of the Northern Pre-Urals[END_REF][START_REF] Calandra | Isotopic partitioning by small mammals in the subnivium[END_REF][START_REF] Eskelinen | Diet of the wood lemming Myopus schisticolor[END_REF].

Contrary to Lemmus, this taxon avoids bogs and marshes and is present throughout the taiga forests from Scandinavia to the eastern coast of the Sea of Okhotsk.

Although Myopus and Lemmus are separated by a large phylogenetic distance [START_REF] Abramson | Genetic analysis of type material of the Amur lemming resolves nomenclature issues and creates challenges for the taxonomy of true lemmings (Lemmus, Rodentia: Cricetidae) in the eastern Palearctic[END_REF] and are easily distinguishable in their external appearance [START_REF] Gromov | Voles (Microtinae). Fauna of the USSR[END_REF]Wilson et al., 2017), they express a great deal of similarity in terms of their dental morphology [START_REF] Borodin | A Diagnostic Guide to Teeth of Arvicolines of the Urals and Western Siberia, from the Late Pleistocene to the Present[END_REF][START_REF] Chaline | Discrimination de la morphologie dentaire de Lemmus lemmus et Myopus schisticolor (Arvicolidae, Rodentia) par l'analyse multivariée[END_REF][START_REF] Koenigswald | Revision of the fossil and recent Lemminae (Rodentia, Mammalia). Papers in Vertebrate Paleontology honoring Robert Warren Wilson[END_REF][START_REF] Kowalski | Lemmings [Mammalia, Rodentia] as indicators of temperature and humidity in the European Quaternary[END_REF], resulting in their fossil remains often being impossible to separate or being lumped together as "Lemmus-Myopus" [START_REF] Kowalski | Pleistocene rodents of Europe[END_REF]. This problem of differentiation has therefore greatly reduced the palaeoecological utility of these taxa, which is particularly problematic considering their occurrence in many Quaternary deposits. Several authors have accordingly attempted to separate these two genera through morphometrical analysis of the first lower molar (m1) and the third upper molar (M3) but with varying degrees of accuracy [START_REF] Chaline | Discrimination de la morphologie dentaire de Lemmus lemmus et Myopus schisticolor (Arvicolidae, Rodentia) par l'analyse multivariée[END_REF][START_REF] Markova | Unification of Criteria for Distinguishing Morphotypes of Cheek Teeth in Lemmings (Lemmini, Arvicolinae, Rodentia)[END_REF]Ponomarev et al., 2013b;[START_REF] Smirnov | Complicated cases of identifying rodent teeth from Late Pleistocene and Holocene deposits of tundra regions of Northern Eurasia[END_REF].

The present paper aims to resolve this problem by providing a robust means of differentiating between the teeth of these two genera, Lemmus and Myopus. This study uses geometric morphometrics, which has never been applied on these two genera, despite its efficiency in identifying similar species in modern and fossil record (e.g. [START_REF] Mcguire | Identifying California Microtus species using geometric morphometrics documents Quaternary geographic range contractions[END_REF][START_REF] Cucchi | The changing pace of insular life: 5000 years of microevolution in the orkney vole (Microtus arvalis orcadensis)[END_REF][START_REF] Hulme-Beaman | Exploring Rattus praetor (Rodentia, Muridae) as a possible species complex using geometric morphometrics on dental morphology[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][START_REF] Navarro | Identifying Past Remains of Morphologically Similar Vole Species Using Molar Shapes[END_REF][START_REF] Kolendrianou | The Palaeolithic cave of Kalamakia (Mani Peninsula), Greece: new insights on the palaeoenvironment using microvertebrates and mesowear analysis of ruminant teeth[END_REF]. Several modern populations were considered, both to establish differences in their dental morphologies and to evaluate the reliability of traditional criteria previously used to separate the two genera. The geometric morphometric approach is then coupled with a discriminant function analysis in order to separate these two genera before being applied to European fossil specimens from two Late Pleistocene sites. The outcomes of this paper consequently propose not only a novel and replicable means of identifying and separating the two genera, but also allow a new view of the past environments and biogeographical evolution of lemming populations.

Taxonomy and Morphology

Lemmus and Myopus genera belong to the tribe Lemmini (Gray, 1825), which is part of the Arvicolinae subfamily (Gray, 1821). The genus Lemmus is currently divided into four distinct species [START_REF] Abramson | Genetic analysis of type material of the Amur lemming resolves nomenclature issues and creates challenges for the taxonomy of true lemmings (Lemmus, Rodentia: Cricetidae) in the eastern Palearctic[END_REF]Wilson et al., 2017): the Nearctic North American brown lemming Lemmus trimucronatus (Richardson, 1825) and the Palearctic Norway lemming Lemmus lemmus (Linnaeus, 1758), Siberian lemming Lemmus sibiricus (Kerr, 1792) and the Amur lemming Lemmus amurensis (Vinogradov, 1924). Originally, Myopus was recognized as a member of the genus Lemmus and named as Lemmus schisticolor (Lilljeborg, 1844) due to osteological similarities, such as the general shape of the skull [START_REF] Ognev | Mammals of the U.S.S.R and adjacent countries. Rodents. Israel Program for Scientific Translation[END_REF], the convergence of the maxillary tooth rows and the inclination of the mandibular tooth rows [START_REF] Abramson | Evolutionary trends in the dentition of true lemmings (Lemmini, Cricetidae, Rodentia): functional-adaptive analysis[END_REF]). Miller (1910, 1912) separated Myopus from Lemmus, leading to the recognition of Myopus schisticolor, a species characterized by a combination of Lemmus characters (skull and teeth) and vole characters (general body-form and foot structure). Recent DNA studies clearly attest to the separation of these two genera [START_REF] Abramson | Genetic analysis of type material of the Amur lemming resolves nomenclature issues and creates challenges for the taxonomy of true lemmings (Lemmus, Rodentia: Cricetidae) in the eastern Palearctic[END_REF][START_REF] Buzan | Mitochondrial phylogeny of Arvicolinae using comprehensive taxonomic sampling yields new insights: Phylogeny of Arvicolinae[END_REF][START_REF] Cook | Historical biogeography at the crossroads of the northern continents: molecular phylogenetics of red-backed voles (Rodentia: Arvicolinae)[END_REF][START_REF] Robovský | Phylogeny of Arvicolinae (Mammalia, Cricetidae): utility of morphological and molecular data sets in a recently radiating clade[END_REF].

Molars of Lemmus and Myopus are hypsodont with a flat grinding surface and have a prismatic structure composed of alternate triangles and deep re-entrant angles filled with cement (fig. 1).

The first lower molar (m1) is constituted by a posterior loop (PL), three closed alternating triangles (T1-T3), two confluent triangles (T4 and T5) and an anterior loop (AL). The third upper molar (M3) shows a similar structure (fig. 1), with an anterior loop, five triangles and a posterior loop (PL). For these two molars, some authors note a high variation in the outer loops (anterior loop for m1, posterior loop for M3), which exhibit different stages of complexity, with the development of supplementary re-entrant angles on both labial and buccal sides (e.g. [START_REF] Markova | Unification of Criteria for Distinguishing Morphotypes of Cheek Teeth in Lemmings (Lemmini, Arvicolinae, Rodentia)[END_REF]Ponomarev et al., 2013b). In general, Lemmus shows a higher tendency to complexity (i.e. development of supplementary re-entrant angles) and Myopus has a simpler morphology. However, due to a large overlap, these trends cannot be used to definitively attribute individuals to one or the other of these genera. Based on biometrical approaches, linear measurements are the most commonly used features to distinguish these two genera on isolated teeth (e.g. [START_REF] Borodin | A Diagnostic Guide to Teeth of Arvicolines of the Urals and Western Siberia, from the Late Pleistocene to the Present[END_REF]Ponomarev et. al., 2013b;Roberts and Parfitt, 1999;[START_REF] Tiunov | The distribution history of the Amur brown lemming (Lemmus amurensis) in the Late Pleistocene -Holocene in the southern Far East of Russia[END_REF], with the following features identified as being of potential significance: 1) differences in the length/width ratio of M3 [START_REF] Smirnov | Complicated cases of identifying rodent teeth from Late Pleistocene and Holocene deposits of tundra regions of Northern Eurasia[END_REF], with the M3 shape in Myopus tending to be more compact and that in Lemmus tending to be more elongated; 2) ratios obtained from several linear measurements to characterize the global shape of the different molars [START_REF] Chaline | Discrimination de la morphologie dentaire de Lemmus lemmus et Myopus schisticolor (Arvicolidae, Rodentia) par l'analyse multivariée[END_REF]; and 3) the size of the teeth, with a greater tooth length in Lemmus [START_REF] Chaline | Discrimination de la morphologie dentaire de Lemmus lemmus et Myopus schisticolor (Arvicolidae, Rodentia) par l'analyse multivariée[END_REF][START_REF] Smirnov | Complicated cases of identifying rodent teeth from Late Pleistocene and Holocene deposits of tundra regions of Northern Eurasia[END_REF]. Finally, a supplementary criterion based on cementum was recently proposed by Ponomarev et al. (2013b). This last work shows that the frequency of cementum deposits in the last re-entrant angles of the outer loops of both M3 and m1 (LRA3 and BRA4 of m1) differs between the two genera. The cementum deposits are absent in Myopus teeth, whereas they are present in 40 to 70% of Lemmus individuals [START_REF] Cheprakov | Встречаемость отложений цемента в терминальных синклиналях моляров у евразийских представителей трибы Lemmini (Rodentia, Arvicolinae)[END_REF][START_REF] Markova | Unification of Criteria for Distinguishing Morphotypes of Cheek Teeth in Lemmings (Lemmini, Arvicolinae, Rodentia)[END_REF]Ponomarev et al., 2013).

Materials and methods

Modern specimens

Since the focus of this paper is the application of the geometric morphometrics method to [START_REF] Abramson | Genetic analysis of type material of the Amur lemming resolves nomenclature issues and creates challenges for the taxonomy of true lemmings (Lemmus, Rodentia: Cricetidae) in the eastern Palearctic[END_REF]. Few individuals of eastern L. sibiricus from Chukotka district were included to capture the range of potential variability of the genus. Additionally, 12 m1 and 8 M3 of North American Brown Lemming (L. trimucronatus) were included to clarify the distinction amongst the different representatives of Lemmus. Numbers and origins of the specimens are detailed in table 1 andfigure 2. Only wellpreserved and complete molars were used in the analyses; all broken and digested teeth were discarded.

Fossil samples

Fossil specimens come from two archaeological/palaeontological sites from Late Pleistocene, la Grotte des Gorges (France) and Gully Cave (England) (fig. 2). As for the modern material, broken and digested teeth were not included in this study. All dates of la Grotte des Gorges have been calibrated using IntCal13 [START_REF] Reimer | IntCal13 and Marine13 radiocarbon age calibration curves 0-50,000 years cal BP[END_REF].

La Grotte des Gorges is located close to Amange village (Jura, France) and was excavated between 2008 and 2017. Two layers (1a and 1b) yielded archaeological material [START_REF] David | La grotte des Gorges ( Jura) : un site inédit à l'interface des territoires symboliques du Paléolithique supérieur ancien[END_REF][START_REF] David | La grotte des Gorges ( Jura) : un site inédit à l'interface des territoires symboliques du Paléolithique supérieur ancien[END_REF], from the end of MIS 3 with AMS dates ranging from 33,030 ± 750 to 29,390 ± 170 uncal. BP. The site yields a diverse association of mammalian species demonstrating a transition from arctic toward steppe environment. Lemmus/Myopus remains come from the upper part of the sequence (layer 1a), dominated by typical taxa of cold-climate phases such as reindeer (Rangifer tarandus), mammoth (Mammuthus primigenius), rhinoceros (Rhinocerotidae) and collared lemmings (Dicrostonyx torquatus) The second phase of the sequence (layer 1b) is dominated by bison (Bison priscus), which is typical taxa of steppe environment. A total of 11 first lower and 3 third upper molars attributed to Lemmus/Myopus has been analysed (table 1 andfigure 2

Landmark schemes and data acquisition

All teeth were photographed using a macro objective (Canon EOS6D mark II, macroscopic between-picture deformation, all specimens were photographed one by one with the same distance from the objective and with the same focus. Outlines were then extracted from pictures using ImageJ v.1.52a. The m1 outline was then orientated with a manual marking of two landmarks (Lm 3 and 7, fig. 1), which define the tooth's vertical axis, in order to correctly position the tooth following the protocol of [START_REF] Navarro | Identifying Past Remains of Morphologically Similar Vole Species Using Molar Shapes[END_REF]. No satisfactory orientation could be obtained with this protocol on the M3, so the procedure was adapted as follows. The outline pixels were rotated along the first principal component axis of the outline pixel coordinate matrix. These standardized orientations then enabled landmarks located at the extreme tips of the salient and reentrant triangles to be automatically detected, as well as those at the extreme tips of the posterior and anterior loops. All landmarked teeth were further inspected for gross error, and the entire procedure was performed by a single operator, thereby avoiding inter-operator biases.

In total, seven landmark schemes and two biometric measurements were defined to describe the morphological variability of the m1 and M3 (table 2 and fig. 1). The global shape of each tooth was analysed with a fixed landmark scheme (FL), using 14 landmarks for the m1 (m1-FL) and 11 landmarks for the M3 (M3-FL). A total of 20 semi-landmarks for m1 and 50 semilandmarks for M3 were used to characterize the outer loops (anterior loop dataset m1-AL and posterior loop dataset M3-PL). The SL schemes, m1-SL and M3-SL correspond to FL datasets combined with semi-landmarks of outer loops (m1-AL or M3-PL). Another semi-landmark scheme describes the alternation between T2 and T3 of M3 (M3-TR) by using a total of 4 fixed landmarks and 30 semi-landmarks. In all semi-landmark configurations, few fixed landmarks were removed to treat whole loops only using semi-landmarks (landmark removed for each dataset: m1-SL/AL: Lm8, 9, 10; M3-SL/PL: Lm 9, M3-TR: Lm3).

In addition, the tooth length (TL) of the two teeth was estimated based on the Euclidean distance between the most anterior and posterior landmarks (Lm1 and 9 for m1; Lm1 and 7 for M3).

Geometric morphometrics

All configurations of each landmark dataset (m1-FL, m1-SL, m1-AL, M3-FL, M3-SL, M3-PL, M3-TR) were superimposed using partial Generalized Procrustes Analysis (GPA), performed independently for each dataset. This procedure aims to separate the shape component from the position, size and orientation components contained within the coordinates of landmark configurations. GPA proceeds by 1) translating each individual configuration so that their centroids coincide, 2) scaling each landmark configuration to unit centroid size, and 3) rotating each configuration onto the dataset mean shape by minimizing Procrustes distances between each configuration and this mean shape [START_REF] Dryden | Statistical shape analysis[END_REF][START_REF] Rohlf | Extensions of the Procrustes Method for the Optimal Superimposition of Landmarks[END_REF][START_REF] Zelditch | Introduction[END_REF]. For m1-SL, m1-AL, M3-SL, M3-PL and M3-TR, semi-landmarks are allowed to slide to minimize the bending energy between the specimen configurations and the meanshape of the dataset under consideration [START_REF] Bookstein | Morphometric Tools for Landmark Data: Geometry and Biology[END_REF].

The obtained coordinates (Procrustes coordinates) were then projected onto the Euclidean tangent space (the tangent point being the mean shape of the dataset). These tangent space coordinates are finally rotated by a Principal Component Analysis (PCA) to keep only principal components with non-null eigenvalues (the rank deficient nature of tangent space data resulting from loss of degree of freedom during GPA, and, in our case, to some datasets with fewer individuals than shape variables). For each dataset, a few outliers were identified as those individuals with high Procrustes distances compared to their mean shape (i.e. lying beyond the third quartile increased by 1.5 times the interquartile range) and were sequentially removed after a new GPA (table 2).

In order to avoid any influence of fossil specimens on the molar shape variation depicted by the tangent space, which was then used to build the discriminant model among both genera, we treated fossil specimens as supplementary specimens in our analyses (see [START_REF] Navarro | Identifying Past Remains of Morphologically Similar Vole Species Using Molar Shapes[END_REF].

Fossil specimens were so aligned on the modern mean shape so that, for a given dataset, they lay on the same shape and tangent spaces as the modern individuals. Then, the rotation computed on the tangent space for the modern individuals was used to rotate the projection of fossil specimens in the same way.

All geometric morphometrics analyses were performed with the geomorph package version 3.2.0 on R v3.6.2 (R Core team, 2019).

Statistical classification

The predictive model was built using linear discriminant analysis (LDA) quantifying the shape differences between the two genera (Myopus vs Lemmus), and the quality of the obtained model was assessed by Leave-One-Out Cross-Validation (LOOCV). This procedure removes one specimen from the sample at a time and predicts its classification between groups using LDA functions calculated on the remaining specimens. Then, the removed specimen is classified by a new function, computed using data that do not take the specimen into consideration, thereby reducing potential inflation of the prediction error and avoiding issues of circular reasoning and over-fitting problems (Kovarovich et al., 2011).

LDA can be extremely sensitive to the relation between the number of groups, the number of predictor variables and the total sample size. As rule-of-thumb, the total number of predictor variables needs to be smaller than the sample size of the smallest predicted group (Kovarovich et al., 2011). Special attention has been paid to this rule due to our small sample size and the high number of predictor variables. Firstly, we reduced the number of predictors used to build the discriminant functions by applying the method of Baylac and Friess (2005). This method substitutes primary data (Procrustes coordinates) by PCs scores to build the model. This was done by building each model and its LOOCV based on the first k PCs, with k varying from 1 to p (p being the sample size of the smallest group under consideration). The optimal number of PCs for the most efficient model was determined by the lowest prediction error rate (i.e. the maximum number of correctly assigned specimens).

Secondly, the method of [START_REF] Evin | The long and winding road: identifying pig domestication through molar size and shape[END_REF] was used to evaluate possible biases due to our unbalanced dataset. For a given unbalanced dataset, 100 balanced samples were built by randomly picking a subset of individuals from the largest group, and for each such sample, LDA was computed. The observed classification score for each unbalanced dataset was then compared with these 100 balanced samples. Models were also compared with the null hypothesis of no shape difference between genera (hence the classification would be only due to chance, resulting in a 50% error rate). 100 such null models were generated for balanced samples of the two species by randomly reassigning specimens to genera. Finally, the quality of the models was evaluated by checking at posterior probabilities of correctly assigned genus for each dataset.

All codes including automatic landmark placement and data treatment are available on request.

Results

Morphospaces and morphological variation

The first lower molar

In the tooth shape analyses (m1-FL, m1-SL and m1-AL), the two first PCs account, respectively, for 32.6%, 56.6% and 69.6% of the explained variance (table 3, fig. 3). The number of PCs explaining at least 90% of the initial shape variation is 13, 12 and 6, respectively.

For the m1-FL configuration, shape changes mainly concern the buccal triangles (T2-T4) shifting labially and the Lm1 and 9 shifting laterally. For the m1-SL, the main behaviour is focused on the anterior loop, which tends to roundness enhancing the confluence of T4/T5 (fig. 3B). The shape variation on m1-AL exhibits a more pronounced development of the re-entrant angles (fig. 3C). The most important part of the variation is expressed by PC1 and tends to separate the two genera for all datasets, with the least overlap for the m1-SL analysis and the largest one for the m1-FL (fig. 3A). On the two first PCs from these three landmark schemes, one can note no specific distinction between Lemmus species (fig. 3). The tooth lengths of Lemmus and Myopus range from 1.6 to 2.4 mm and from 1.4 to 1.7 mm, respectively (fig. 3H).

A large overlap in the tooth length range is observed between the two genera, with 86% of the Myopus specimens being included within Lemmus's tooth length range. Within the Lemmus genus, L. sibiricus is the species which has the longest tooth length compared to the two others.

The third upper molar

The same pattern as m1 is observed on the M3 with the first two PCs explaining 46.1%, 70.0%, 75.7% and 57.7% of variation respectively M3-FL, M3-SL, M3-PL and M3-TR (table 3, fig.

3). The number of PCs explaining at least 90% of the initial shape variation is 10, 8, 5 and 6, respectively. The shape variation illustrated by PC1 for M3-FL concerns the width and length of the teeth, with a tendency to be more "compact" (i.e. thinner and smaller) (fig. 3D). For M3-SL, the main changes occur on the lingual triangles and anterior loop, with an enlargement on the lingual side (fig. 3E). Similar changes on the posterior loop are also observed on M3-PL.

The first PC of M3-TR describes the T3/T4 alternation, with a T3 tending to close. As for m1, the PC1 tends to separate the two genera, with least overlap for the M3-SL analysis and the largest overlap for the M3-TR for which no separation is observed between the two genera (fig. 3G). Again, the two first PCs do not show any clear difference between Lemmus species across all landmark schemes. The M3 tooth lengths follow the same pattern as for m1 with a large overlap between the two genera. Lemmus range from 1.2 to 1.8 mm and Myopus from 1.0 to 1.4 mm, with 90% of M. schisticolor being included in the Lemmus tooth length interval (fig. 3I).

Statistical classification of modern specimens and shape differences

Although there was only a relatively small number of individuals from Myopus (the smallest sample), the unbalanced prediction error curves of the landmark schemes do not show particular bias compared to the balanced ones (fig. 4a). Nonetheless, a bias was noted between the dimensionality of the predictor (the number of PCs included to build the LDA) and the prediction error, with an increase when a high number of predictors is used (fig. 4a). This pattern could be due to the small size of the model samples in relation to the large dimensionality of the shape space, as encountered by [START_REF] Evin | The long and winding road: identifying pig domestication through molar size and shape[END_REF]. It is likely that with a larger dataset, such as that presented by [START_REF] Navarro | Identifying Past Remains of Morphologically Similar Vole Species Using Molar Shapes[END_REF], this pattern would not be observed. The optimal number of predictors (PCs), i.e the one giving the lowest prediction error rate, was selected to build the discriminant functions of the models and are detailed in table 3. Six models (m1-FL, m1-SL, m1-AL, M3-FL, M3-SL and M3-PL) have low prediction errors inferior to 5% for unbalanced samples (the real samples). Models of the m1 tend to have a lower prediction error compared to those from the M3 (fig. 4b; table 3). Landmark schemes of the m1 taking in account the anterior loop (m1-SL and m1-AL) give the best classifications, up to 0.96 and 0.93. The same pattern is visible on the M3 and the posterior loop, with M3-SL having 0.95 optimal classification score, with M3-FL and M3-PL having comparable scores of 0.89 and 0.88. In contrast, the tooth length (m1-TL and M3-TL) has very low classification scores, with 0.41 and 0.42 respectively, which is comparable to a classification only due to chance.

Figure 5 illustrates the shape differences between Lemmus and Myopus for modern individuals along the first discriminant axis from the LDA regarding the different landmark schemes. In terms of the m1 shape, Myopus tends to differ from Lemmus with: 1) a different orientation of the lingual triangles, tending to sweep upwards; 2) a simpler anterior loop with a less marked labial peak; and 3) a LRA3 that extends less deeply over the BRA2 without generating a new re-entrant angle on the buccal side of the anterior loop. For the M3, Myopus tends to differ from Lemmus as follows: 1) an overall stockier tooth; 2) a simpler posterior loop, almost flat; and 3) a shorter distance between Lm2 and Lm8. For m1 and M3, tooth lengths show a large overlap with almost all Myopus plotting with the smaller Lemmus specimens.

Fossil identifications

The best performing model (m1-SL) obtained with modern specimens allows the identification of at least one Myopus individual in Grotte des Gorges and one in Gully Cave from the Allerød interstadial-Younger-Dryas transition level, with posterior probabilities at 0.998 and 0.935, respectively (table 4, figure 6a). For these two individuals, similar attributions to Myopus were obtained with other models (m1-FL, m1-AL) (see supplementary table A). The m1-TL model (table 4, figure 6b) cannot attribute any individuals as a Myopus with significant posterior probabilities and only 9 individuals are confidently identified as Lemmus. Most of the specimens are undetermined.

For the M3 (table 4, figure 6c), one tooth has been attributed to Myopus at Gully Cave for the same level, with the best model (M3-SL) giving 0.994 of posterior probability. This attribution was equally obtained with M3-FL and M3-PL models. No model suggests the presence of Myopus in the M3 sample from Grotte des Gorges, but only three teeth were analysed. As for m1-TL, the M3-TL model cannot attribute any individual to Myopus with high posterior probabilities (figure 6d), and the model is unable to confidently identify a large majority of the specimens (table 4, figure 6).

Discussion

Myopus schisticolor and Lemmus sp. are two lemmings that are currently widespread throughout north Eurasia. They have two distinct habitat preferences, the first inhabits exclusively taiga forest, whereas the second prefers an open tundra environment (e.g. [START_REF] Mitchell-Jones | The atlas of European mammals[END_REF][START_REF] Niethammer | Handbuch der Säugetiere Europas. Band 2/I. Nagetiere II[END_REF]. They both consume mosses, but with variable seasonal proportions: the diet of Myopus schisticolor is restricted only to this vegetation type, whereas Lemmus sp. diversifies its diet during spring and summer with dicots, grasses and sedges being consumed (e.g [START_REF] Batzli | Nutritional Ecology of Microtine Rodents: Food Habits of Lemmings near Barrow, Alaska[END_REF][START_REF] Calandra | Isotopic partitioning by small mammals in the subnivium[END_REF][START_REF] Eskelinen | Diet of the wood lemming Myopus schisticolor[END_REF]Rodger and Lewis, 1986). If distinguishing between these two species is relatively easy at the present day, due to their distinctive external morphology, they are extremely difficult to identify with fossil material, on account of their similar tooth morphologies.

Molar shape differentiation

Despite the similarity in tooth morphology of these two lemming genera, the geometric morphometric approach used in this paper has been able to describe and quantify intergeneric shape differentiation of the first lower molar (m1) and the third upper molar (M3) of modern specimens. Describing the whole shape of the molars with a set of landmarks associated with semi-landmarks on the anterior loop of m1 or posterior loop of M3 (m1-SL and M3-SL) appears to be the best landmark scheme among those tested here for the modern reference framework. These two landmark schemes have allowed us to separate the two genera with high statistical probabilities, around 0.95 confidence (table 3).

These landmark schemes suggest that the main differences between these two genera relate to the tooth complexity with the outer loops more developed in Lemmus than in Myopus. The anterior loop on the m1 or the posterior loop on the M3 tend to exhibit a supplementary angle, which is not observed in Myopus. For Myopus, teeth tend to be stockier and simpler in shape.

This increased complexity of the anterior m1 loop in Lemmus, always absent in Myopus, and the organisation of the buccal triangles shifting anteriorly on Lemmus, are congruent with the taxonomic literature [START_REF] Chaline | Discrimination de la morphologie dentaire de Lemmus lemmus et Myopus schisticolor (Arvicolidae, Rodentia) par l'analyse multivariée[END_REF][START_REF] Markova | Unification of Criteria for Distinguishing Morphotypes of Cheek Teeth in Lemmings (Lemmini, Arvicolinae, Rodentia)[END_REF][START_REF] Smirnov | Complicated cases of identifying rodent teeth from Late Pleistocene and Holocene deposits of tundra regions of Northern Eurasia[END_REF]Ponomarev et al., 2013b). Some of these authors also noted that on the M3, the T3/T4 alternation and asymmetry could be used to distinguish both genera. The morphometrical approach developed in this paper thus offered the opportunity to evaluate these criteria (M3-TR). Despite the relevance of these criteria, we demonstrate here that they cannot be used to clearly differentiate the two genera because of a large overlap in their morphology (figure 3; table 3).

Among the external morphological criteria allowing differentiation between these two genera, their differences in sizes have been regularly mentioned, with Myopus being smaller [START_REF] Borodin | A Diagnostic Guide to Teeth of Arvicolines of the Urals and Western Siberia, from the Late Pleistocene to the Present[END_REF]. This criterion was often applied to fossil material through measurement of the length of the teeth or dental row (e.g. [START_REF] Borodin | A Diagnostic Guide to Teeth of Arvicolines of the Urals and Western Siberia, from the Late Pleistocene to the Present[END_REF][START_REF] Kowalski | Fossil lemmings [Mammalia, Rodentia] from the Pliocene and early Pleistocene of Poland[END_REF][START_REF] Rhodes | Small mammal taxonomy, taphonomy, and the paleoenvironmental record during the Middle and Upper Paleolithic at Geißenklösterle Cave (Ach Valley, southwestern Germany)[END_REF]. Nonetheless, as illustrated in figure 3, there is a high overlap between the two genera in terms of tooth length. This overlap is mainly observed between M. schisticolor and L. lemmus, whereas L. sibiricus teeth are taller, with a clear separation from other lemmings (figure 3H,I). Applying the criterion of size to Western European fossil assemblages with lemmings is therefore likely to conflate L. lemmus and M. schisticolor, affirming presence of the former but obscuring the presence of the latter. In any case, size always has to be used very cautiously on fossil species because it can be highly variable according to the time period, the climate or even the latitude and altitude (e.g. [START_REF] Delpech | Biomasse d'Ongulés au Paléolithique et inférences sur la démographie[END_REF][START_REF] Klein | Carnivore size and quaternary climatic change in southern Africa[END_REF][START_REF] Klein | Glacial/interglacial size variation in fossil spotted hyenas (Crocuta crocuta) from Britain[END_REF].

By its ability to capture the whole geometry of an object, geometrics morphometrics has therefore become increasingly popular in archaeological and palaeontological studies (e.g. [START_REF] Cucchi | The changing pace of insular life: 5000 years of microevolution in the orkney vole (Microtus arvalis orcadensis)[END_REF][START_REF] Escudé | Reappraisal of 'chronospecies' and the use of Arvicola (Rodentia, Mammalia) for biochronology[END_REF][START_REF] Evin | The long and winding road: identifying pig domestication through molar size and shape[END_REF][START_REF] Killick | Geometric Morphometric analysis of the Microtus M1 and its application to Early Middle Pleistocene in the UK[END_REF][START_REF] Marr | Faunal response to abrupt climate change: the history of the British mammal fauna from the Lateglacial to the Holocene[END_REF][START_REF] Mcguire | Identifying California Microtus species using geometric morphometrics documents Quaternary geographic range contractions[END_REF][START_REF] Miele | Deep learning for species identification of modern and fossil rodent molars[END_REF][START_REF] Navarro | Identifying Past Remains of Morphologically Similar Vole Species Using Molar Shapes[END_REF]. Associated with classification analysis, such as a LDA, it is a powerful approach that allows differentiation of morphologically close species, thereby assisting in the robust identification of fossil remains and in our understanding of changing faunal communities.

Myopus fossil identifications: palaeobiogeographical implications

During Late Pleistocene glaciations, cold environmental conditions were favourable to lemmings and facilitated their dispersal throughout a large part of Europe [START_REF] Markova | Institute of geography of Russian Academy of Sciences[END_REF][START_REF] Royer | Etude paléoenvironnementale et paléoclimatique du Pléistocène supérieur du Sud-Ouest de la France, à partir d'analyses comparées d'associations fauniques et de biogéochimies effectuées sur les micromammifères[END_REF], from the Urals (Ponomarev et. al., 2013) to Ireland [START_REF] Monaghan | Irish Quaternary Vertebrates[END_REF][START_REF] Sutcliffe | Pleistocene rodents of the british isles[END_REF]Woodman et al., 1996), including for example France (Marquet, 1993;[START_REF] Royer | Investigating the Influence of Climate Changes on Rodent Communities at a Regional-Scale (MIS 1-3, Southwestern France)[END_REF], Denmark [START_REF] Bennike | A multi-disciplinary macrofossil study of Middle Weichselian sediments at Kobbelgård, Møn, Denmark[END_REF]Larsen and Mangerud, 1990) or Czechoslovakia (Horacek and Sanchez-Marco, 1984). Until the present paper, the presence of Myopus in Western Europe was either not envisaged at all or, at best, still questioned. Our results attest its presence in two sites from two distinct periods of Late Pleistocene: la Grotte des Gorges, located in East of France, which is from the end of MIS 3, around 30 ka cal BP, and Gully Cave, situated in the Somerset (England), from the Allerød interstadial -Younger Dryas transition level.

Occurrences of M. schisticolor have already been suggested for European Early and Middle

Pleistocene site at Boxgrove (England, Roberts and Parfitt, 1999) and Vergranne (France, [START_REF] Chaline | Les lemmings fossiles du genre Lemmus (Arvicolidae, Rodentia)dans le pléistocène de France[END_REF] and suspicions of its presence (records of Myopus/Lemmus) have raised in Chlum-4 (Czech Republic, Kowalski, 2001), Nyaravai-2 (Lithuania, [START_REF] Kowalski | Pleistocene rodents of Europe[END_REF], and Sudmer-Berg 2 (Germany, von [START_REF] Koenigswald | Sudmer-Berg-2, a fauna of the early middle Pleistocene from the Harz[END_REF]. In general, these suspicions are based on both the morphology of the teeth and on the faunal associations, which are characteristic of a forest environment. However, these tentative suggestions were not uniformly accepted (van [START_REF] Kolfschoten | On the application of fossil mammals to the reconstruction of the palaeoenvironment of northwestern Europe[END_REF]1996;[START_REF] Kowalski | Lemmings [Mammalia, Rodentia] as indicators of temperature and humidity in the European Quaternary[END_REF][START_REF] Kowalski | Pleistocene rodents of Europe[END_REF]; for example at Vergranne, [START_REF] Kowalski | Lemmings [Mammalia, Rodentia] as indicators of temperature and humidity in the European Quaternary[END_REF] stated that the Myopus identification cannot be accepted as definitive. His opinion was that it is not possible to separate Lemmus and Myopus in the fossil record, and that Myopus did not reach Europe during Pleistocene [START_REF] Kowalski | Pleistocene rodents of Europe[END_REF]. Approaches as geometric morphometrics or ancient DNA, which provide reliable identification, imply reconsidering certain determinations based on presuppositions of past studies (geographical area, morphological features, faunal assemblages). The same issues were recently underlined for other species as for example Sorex, with the identification of the boreal Sorex tundrensis in Germany during the end of MIS3 [START_REF] Freund | Sesselfelsgrotte I: Grabungsverlauf und Stratigraphie[END_REF][START_REF] Prost | Effects of late quaternary climate change on Palearctic shrews[END_REF], a species that was not supposed to reach Europe.

Lemming specimens from older Early and Middle Pleistocene sites are generally assigned to ancestral forms of Lemmus, possibly Lemmus kowalskii. This extinct species is then assumed to have lived under temperate conditions in a more wooded environment, as suggested by the forest taxa with which it was associated (van [START_REF] Kolfschoten | On the application of fossil mammals to the reconstruction of the palaeoenvironment of northwestern Europe[END_REF]1996;[START_REF] Kowalski | Lemmings [Mammalia, Rodentia] as indicators of temperature and humidity in the European Quaternary[END_REF]2001). However, identification of Lemmus kowalskii is highly controversial, since it is morphologically similar to Myopus, but much closer in size to Lemmus [START_REF] Harrison | Occurrence of Lemmus kowalskii Carls and Rabeder, 1988 (Rodentia: Microtinae: Lemmus) in the Lower Pleistocene of East Anglia[END_REF].

For the Late Pleistocene, several studies attest the presence of Myopus in the Ural (Ponomarev et al., 2013a, b, c) by using the ratios of [START_REF] Smirnov | Complicated cases of identifying rodent teeth from Late Pleistocene and Holocene deposits of tundra regions of Northern Eurasia[END_REF] or morphotype differentiation associated with non-metric multidimensional scaling (Ponomarev et al., 2013a, c). In Europe, Lemmus / Myopus remains from this period are always identified as Lemmus (Lemmus lemmus in the western part of Europe and Lemmus sibiricus in the eastern part). To our knowledge, only one previous occurrence of Myopus schisticolor has been suggested for this period in western part of Europe, in Level 4 of King Arthur's Cave (Wye Valley, England) [START_REF] Price | Late Pleistocene and Early Holocene Small Mammals in South West Britain: Environmental and Taphonomic Implications and Their Role in Archaeological Research[END_REF], which has been attributed to the Bølling-Allerød interstadial (Bronk Ramsey, 2002). Nevertheless, these determinations (obtained on 12 teeth from a total sample of 17) were only based on tooth length [START_REF] Price | Late Pleistocene and Early Holocene Small Mammals in South West Britain: Environmental and Taphonomic Implications and Their Role in Archaeological Research[END_REF] and were not corroborated by another approach, although the presence of boreal species in the same level, such as Clethrionomys rufocanus, tends to strengthen this hypothesis [START_REF] Price | Late Pleistocene and Early Holocene Small Mammals in South West Britain: Environmental and Taphonomic Implications and Their Role in Archaeological Research[END_REF]. The identification of Myopus in the Allerød -interstadial layer of Gully Cave corroborates Price's suppositions, suggesting the presence of this species during milder climatic phases in the southern part of England. This presence throughout the septentrional region of Western Europe (i.e. latitude > 45°N) during the Late Pleistocene, raises the question of whether its occurrence was (semi)continuous, or sporadic related to specific climatic events favouring its expansion from eastern regions.

Palaeoenvironmental implications of Myopus occurrence in Late Pleistocene

Since Myopus is today closely linked to the taiga, its identification has strong implications for palaeoenvironmental reconstructions, suggesting past local environments with boreal forest cover mainly composed of Pinus, Betula and Picea, which is favourable to the development of the thick moss cover that constitutes nearly all of its diet [START_REF] Bobretsov | Population dynamics of wood lemming (Myopus schisticolor) in different landscapes of the Northern Pre-Urals[END_REF][START_REF] Eskelinen | Diet of the wood lemming Myopus schisticolor[END_REF]. Our results from both sites (Grotte des Gorges and Gully Cave) suggest the association of individuals from both genera Lemmus and Myopus. While this might appear contradictory due to their distinct environmental preferences, they were found within a diversified faunal association, including for instance Dicrostonyx torquatus, Alexandromys oeconomus or Lasiopodomys gregalis. Although these fossil faunal communities are 'disharmonious' by comparison to present-day small mammal associations and might therefore reflect a palimpsest of distinct phases, the stratigaphical integrity of our study samples is clear.

It is therefore more likely that these diverse assemblages reflect the mosaic environment of boreal forest and steppes that are increasingly recognised as typical of the Late Pleistocene of northern Europe.

Multi-proxy investigations on sediment archives have demonstrated continental vegetation to be highly sensitive to Dansgaard-Oeschger events [START_REF] Fletcher | Millennial-scale variability during the last glacial in vegetation records from Europe[END_REF], consequently impacting small mammal communities. Pollen from the Bergsee lacustrine record in Germany demonstrates a high frequency succession of steppe and boreal forest phases, consistent with stadial-interstadial oscillations, and underpinned by increases in Juniperus and Pinus pollen during Greenland Interstadials 8 and 7 [START_REF] Becker | Paleoenvironmental studies on lake Bergsee, Black Forest, Germany[END_REF][START_REF] Duprat-Oualid | Vegetation response to abrupt climate changes in Western Europe from 45 to 14.7k cal a BP: the Bergsee lacustrine record (Black Forest, Germany)[END_REF]. The presence of Myopus in la Grotte des Gorges could therefore be related to one of these Greenland Interstadial events.

For Gully Cave, two levels have been investigated, one attributed to the Allerød interstadial -Younger Dryas transition, and one from the Younger Dryas. Among individuals from the older level, at least two teeth suggest the presence of M. schisticolor according to the most efficient models (m1-SL: 1, M3-SL: 1), but for the cold-climate level of Younger Dryas, no specimens have been attributed to M. schisticolor. Vegetation in Britain was highly sensitive through Late Glacial stadials and interstadials, with regional variation depending of the latitude (e.g. [START_REF] Birks | To what extent did changes in July temperature influence Lateglacial vegetation patterns in NW Europe?[END_REF][START_REF] Huntley | An atlas of past and present pollen maps for Europe, 0-13,000 years ago[END_REF][START_REF] Jones | Pleistocene Environments in the British Isles[END_REF][START_REF] Pennington | The Late Devensian flora and vegetation of Britain[END_REF][START_REF] Walker | Devensian Lateglacial environmental changes in Britain: a multiproxy environmental record from Llanilid, South Wales, UK[END_REF]. In southern England, the Allerød interstadial was a temperate-climate phase characterized by an increase in forest taxa (Hill et al., 2007;[START_REF] Walker | Devensian Lateglacial environmental changes in Britain: a multiproxy environmental record from Llanilid, South Wales, UK[END_REF], in particular the spread of Betula [START_REF] Birks | To what extent did changes in July temperature influence Lateglacial vegetation patterns in NW Europe?[END_REF]. The development of this semi-open boreal environment in southwestern England (Hills et al., 2007) is echoed by the apparent disappearance of Equus ferus at this time [START_REF] Kaagan | The horse in late Pleistocene and Holocene Britain[END_REF], and favoured the coexistence of both arctic and boreal species, such as Cervus elaphus and Rangifer tarandus in Cheddar Gorge [START_REF] Currant | The Mammal Faunas of the British Late Pleistocene[END_REF]. King Arthur's Cave, which is located 80 km north of Gully Cave, has equally yielded both forested and cold-adapted taxa (e.g. Dicrostonyx torquatus, Alexandromys oeconomus, Lasiopodomys gregalis, Clethrionomys rufocanus, Microtus agrestis, Lemmus lemmus; [START_REF] Price | Late Pleistocene and Early Holocene Small Mammals in South West Britain: Environmental and Taphonomic Implications and Their Role in Archaeological Research[END_REF], suggesting a mosaic environment. This environment was able to support both Lemmus and Myopus populations during the Allerød interstadial, with the presence of both open and forested areas. The abruptness of the succeeding Younger Dryas cooling event led to drastic vegetation changes, with a decline in forested environments, a decrease in Betula, and a rise of Artemisia in southern England [START_REF] Birks | To what extent did changes in July temperature influence Lateglacial vegetation patterns in NW Europe?[END_REF].

European Betula macrofossils from this period only belong to Betula nana [START_REF] Birks | To what extent did changes in July temperature influence Lateglacial vegetation patterns in NW Europe?[END_REF], reflecting the reduction of forest and inhibiting the development of the moss cover needed by Myopus.

Conclusion

Our results highlight that: i) the presence of Myopus schisticolor is confirmed for the first time in the Late Pleistocene fossil record of western Europe, although it remains to be established whether its occurrence was (semi)continuous or sporadic, only related to specific climatic events favouring its expansion from eastern regions; ii) its presence has important consequences for palaeoenvironmental interpretations, implying the existence of boreal open/semi open environments. Muonio,4. Rovaniemi,5. Sabetta,6. Erkuta,7. Bely Island,10. Sotkamo area,11. Posio,12. Grotte des Gorges, 13. Gully Cave. Localities from North America and Eastern Siberia are not shown (see text for detail).

Figure captions

  ).GullyCave is located in Ebbor Gorge (Somerset, UK) and has been the subject of ongoing excavation since 2006. The upper stratigraphic units (of relevance here) consist of matrix-rich limestone breccias of Lateglacial Interstadial and Holocene age respectively, separated by a coarse, open-framework limestone breccia deposited during the Younger Dryas, the stratigraphy underpinned by radiocarbon dates on mammalian faunal remains from throughout the sequence. The cave was completely filled by around 10,000 cal years BP and the sequence was capped by a discontinuous flowstone. A total of 35 m1 teeth were analyzed, 12 from the Younger Dryas levels (date range 12409-12037 cal. BP to 11950-11259 cal. BP) and the other 23 from the boundary level between the Allerød interstadial and the Younger Dryas (date range 12379-12637 cal. BP to 13214-13444 cal. BP), and 43 M3s, with 17 from the Younger Dryas and 26 from the Allerød interstadial-Younger Dryas transition. The mammalian assemblage from the Younger Dryas levels sampled here includes Lepus timidus, Myodes sp., Dicrostonyx torquatus, Lemmus/Myopus, Microtus agrestis or arvalis, Lasiopodomys gregalis, Alexandromys oeconomus, Vulpes sp. and Rangifer tarandus. The sample from the Allerød interstadial-Younger Dryas transition has yielded L. timidus, Ochotona pusilla, Arvicola amphibius, Clethrionomys glareolus, D. torquatus, Lemmus/Myopus, M. agrestis or arvalis, L. gregalis, Al. oeconomus, V. lagopus, R. tarandus and Cervus elaphus.

Figure 1 :

 1 Figure 1: Landmark schemes (in grey squares) defined for Lemmini first lower molar (m1) and third upper molar (M3). Hand-placed landmarks are in white, automatically detected landmarks are in black, and semi-landmarks are in blue (for m1-SL, m1-AL, M3-SL and M3-AL) and in yellow (for M3-TR). Abbreviations: T = triangles, LRA = lingual re-entrant angles, BRA = buccal re-entrant angles, AL = anterior loop, PL = posterior loop. Modified after van der Meulen (1973).

Figure 2 :

 2 Figure 2: A) Global climate evolution and chronology of the last part of the last glacial following the NGRIP-ice record (Andersen et al., 2004) with the position of the two fossil sites. Abbreviations: MIS: marine isotope stage, HE: Heinrich events, LGM: last glacial maximum, BA: Bølling -Allerød interstadial, YD: Younger-Dryas. B) Geographical range of Lemmus lemmus (blue), Lemmus sibiricus (red) and Myopus schisticolor (yellow) with the overlap area (orange). Circles indicate the location of the modern samples and squares (brown) the localisation of the fossil samples. Site number: 1. Karigasniemi, 2. Kilpisjärvi, 3. Pällasjarvi-

Figure 3 :

 3 Figure 3: Morphospaces (PC1 -PC2 planes) of modern specimens with associated shape variation for all landmark schemes (A: m1-FL; B: m1-SL; C: m1-AL; D: M3-FL; E: M3-SL;F: M3-PL) and violin plots for tooth length (H: m1-TL; I: M3-TL). For PCA, displayed shape variation corresponds to morphological changes along the 1st PC axis, depicted by lollipop graphs (grey outline and black dots for the mean shape, black segments for the shape deviation from the mean shape to the most distant individual score on the considered PC). Percentage of shape variation explained by each PC are reported. Myopus are in yellow, Lemmus lemmus in blue, Lemmus sibiricus in red and Lemmus trimucronatus in cyan.

Figure 4 :

 4 Figure4: A) Variation of the LDA's prediction error according to the number of PCs included to build the model for each landmark scheme, with m1 on the left and M3 on the right. For each plot, continuous coloured lines show the prediction error for the original samples, dotted coloured lines show the prediction error for the set of 100 balanced samples, and the grey line show the prediction error for the set of 100 randomly reassigned samples with the associated error bars standing for the 95% confidence intervals (null hypothesis). B) Model prediction error with the optimal number of PCs for each landmark scheme. Coloured points correspond to the observed prediction error on the real datasets, black points correspond to the prediction error for the set of 100 balanced samples with the associated 95% confidence interval, grey points show the prediction error for the set of 100 randomly reassigned samples with the associated error bars standing for the 95% confidence interval (null hypothesis).

Figure 5 :

 5 Figure 5: Distribution histograms of the modern samples along the linear discriminant axis with associated shape changes. Lemmus are in blue, Myopus are in yellow. The displayed shapes correspond to the extreme individual scores on the discriminant axis for each landmark scheme.

Figure 6 :

 6 Figure 6: Distribution histograms of the fossil samples along the discriminant axis from the model based on modern samples for m1-SL, m1-TL, M3-SL and M3-TL datasets. Fossil samples are reported by archaeological/palaeontological levels (scatter plot). Lemmus are in blue, Myopus are in yellow.

  FL = Fixed landmarks, SL= Semi landmarks (including FL and AL for m1 and FL and PL for M3), PL = Posterior loop, AL = Anterior loop, TL = ToothLength. Table 3: Abbreviations: Tot. PCSs= Total Number of PCs, Opt. PCs = Optimal Number of PCs used in the model, PC1+PC2 = cumulative proportion of variance explained by PC1 + PC2, % Opt. class. actual = percentage of prediction error of modern specimens with the optimal number of PCs.

  western European fossil sites, modern Lemmus and Myopus individuals were targeted from Scandinavia and Russia. The modern samples are composed of research laboratory collections and Snowy Owl (Bubo scandiacus) pellets and are housed in the research laboratory collections in Biogéosciences at the University of Burgundy and in the Arctic Research Station of Labytnangi in Russia. In total, 30 Myopus and 96 Lemmus m1 specimens, and 31 Myopus and 91 Lemmus M3 specimens were analysed (see table 1 for details). Since the Myopus genus is monospecific, our sample was only composed of individuals of Myopus schisticolor. To describe most of the current morphological variability of the Lemmus genus, our sample included specimens from the three different species. Most of our sample belong to two

	palearctic species, with 50 Norway Lemming (L. lemmus) and 34 western Siberian Lemming
	(L. sibiricus). The two species are closely related in their phylogeny but were already separated
	during the Late Pleistocene

Table captions Table 1 :

 captions1 List of the analysed specimens with their origin and counts (neighbour localities are merged). See figure1for localities.

Table 2 :

 2 Datasets used in the present study.

Table 4 :

 4 Fossil identification for m1 and M3. The number of specimens identified with associated posterior probabilities > 0.9 are in brackets

Table 1 :

 1 List of the analysed specimens with their origin and counts (neighbouring localities are merged). See figure1for localities.

	948
		Species	Location	Lower m1	Upper M3
	Modern specimens	Lemmus lemmus (50)	1-Karigasniemi	3
			2-Kilpisjärvi	15
			3-Pallasjärvi, Muonio	11
			4-Rovaniemi	21
		Lemmus sibiricus (34)	5-Sabetta	9
			6-Erkuta	10
			7-Bely Island	11
			(8)-Wrangel Island	4
		Lemmus trimucronatus (12)	(9)-Victoria Island	12
		Myopus schisticolor (31)	10-Sotkamo, Valtimo, Kuhmo, Lieka,	20

Table 2 :

 2 Datasets used in the present study. FL = Fixed landmarks, SL= Semi landmarks (including FL and AL for m1 and FL and PL for M3), PL = Posterior loop, AL = Anterior loop, TL = ToothLength.

		Landmark scheme	Number of Landmarks Semi-landmarks	Evaluated criteria	Reference
		FL	11	-	Width of the teeth	Chaline, 1988
		PL	2	50	Complexity of the posterior loop	Ponomarev et al., 2013, Markova et al., 2017
		SL	9	50	Complete shape	-
	M3				Depth of A1	Chaline, 1988
		TR	4	10 + 20	Confluence of FL2	Ponomarev et al., 2013
					Asymmetry of FL2	Markova et al., 2017
		TL	2	-	Tooth length	Chaline, 1988
		FL	14	-	Organisation of the buccal triangles Chaline, 1988
		AL	2	20	Complexity of the anterior loop	Ponomarev et al., 2013, Markova et al., 2017
	M1	SL	11	20	Complete shape	-
		TL	2	-	Tooth length	Chaline, 1988

Table 3 :

 3 Abbreviations: Tot. PCSs= Total Number of PCs, Opt. PCs = Optimal Number of PCs used in the model, PC1+PC2 = cumulative proportion of variance explained by PC1 + PC2, % Opt. class. actual = percentage of prediction error of modern specimens with the optimal number of PCs (see text for details).

									Fossil
		Dataset	Tot. PCs	Opt. PCs	PC1 + PC2	% Pred. Error	% Opt. class. actual	% Opt. class. fossil	identification Lemmus Myopus
		FL	28	14	32.6	2.5	85.5	80.4	38	8
		SL	62	19	56.6	0.8	95.9	97.8	44	2
	M1	AL	44	12	69.6	3.3	93.4	71.7	27	19
		TL	-	-	-	19.8	41.3	19.6	39	7
		FL	18	17	46.1	4.1	89.3	80.4	36	1
		SL	99	26	70	3.4	94.9	100	42	1
	M3	PL	92	3	75.7	3.4	88.1	65.2	22	1
		TR	62	19	57.7	7.3	83.6	76.1	27	5
		TL	-	-	-	20.5	41.8	34.8	13

Acknowledgments

This study was partially supported by the project HARCGLOB (AAP 2020 Région Bourgogne Franche-Comté). The authors thank Olivier Gilg for giving access to specimens from Victoria Island. We also thank Jérôme Thomas and ReColNat for giving access to the lemming specimens housed at Biogeosciences Palaeontology collection, as well as Aleksandr Sokolov and Natalia Sokolova of the Ural branch of Russian Science Academy for their loan of Siberian lemming's specimens. Danielle Schreve acknowledges the support of the National Trust and Natural England in the excavation of Gully Cave.

Author contributions

Conceptualization, L.A., S.M. and A.R.; Methodology and Software, L.A. and R.L.; Formal Gully Cave (Younger Dryas)

Lemmus 12 ( 12) 12 ( 12) 9 ( 2) 9 ( 0) 15 ( 15) 17 ( 16) 15 ( 10) 16 ( 1) 14 ( 2)

Gully Cave (Allerød-Younger Dryas)

Myopus 6 (4) 1 ( 1) 13 ( 9) 3 ( 0) 2 ( 1) 1 ( 1) 10 ( 6) 3 ( 1)

Lemmus 17 ( 10) 22 ( 22) 10 ( 9) 20 ( 3) 24 ( 19) 25 ( 25) 16 ( 9) 23 ( 19) 25 ( 12)

Grotte des Gorges Myopus 2 (2) 1 ( 1) 3 (3) 1 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)

Lemmus 9 (9) 10 ( 9) 8 ( 7) 10 ( 6) 3 (2) 3 (3) 3 (3) 3 (3) 3 (3)


Simon Ladouce

Juan Jesus

Torre Tresols

Ludovic Darmet

Frédéric Giuseppe Ferraro

Improving user experience of SSVEP-BCI through reduction of stimuli amplitude depth

Improving user experience of SSVEP-BCI through reduction of stimuli amplitude depth S. Ladouce 1 , J. J. Torre Tresols 1 , L. Darmet 1 , G. Ferraro 1 , and F. Dehais 1,2 Abstract-Steady-States Visually Evoked Potentials (SSVEP) is currently one of the most widely used paradigms in Brain-Computer Interfaces (BCI). Although the SSVEP-BCI are characterized by their high and robust classification performance, the repetitive presentation of flickering stimuli is uncomfortable from a user experience perspective point of view. Indeed, the low-level visual features of SSVEP stimuli make them straining to the eyes over time and could be disruptive to the execution of tasks requiring sustained attention. They could even induce epileptic seizures. This study explores the reduction of the stimulation amplitude depth (a magnitude diminution of 90%) for the design of SSVEP stimuli as a solution to improve user comfort. The classification accuracy obtained by different pipelines was systematically compared between low amplitude and standard full amplitude SSVEP stimuli. The results reveal high classification accuracy for both high (99.8%) and low magnitude (80.2%) stimuli using the Task-Related Component Analysis (TRCA) classification method. The present findings demonstrate the validity of reducing SSVEP stimuli amplitude to increase users' comfort paving the way for transparent BCI operation.

I. INTRODUCTION

Steady-States Visually Evoked Potentials (SSVEP) refer to the sustained cortical entrainment to the frequency of periodical visual stimulations. When attending to repetitive visual stimuli (RVS) whose low-level features (color, luminance, shape) fluctuate at a certain rate (commonly through flickering shapes and alternating checkerboards), neuronal populations in the visual cortex synchronize their firing rate to the frequency of the visual stimulation [START_REF] Vialatte | Steady-state visually evoked potentials: Focus on essential paradigms and future perspectives[END_REF]. The SSVEP is therefore well localized in the frequency domain of the electroencephalography (EEG) signal. Its signal in the EEG is most prominent at occipital sensor sites and source locations [START_REF] Srinivasan | Steady-state visual evoked potentials: Distributed local sources and wave-like dynamics are sensitive to flicker frequency[END_REF]. The SSVEP has been reported over a wide range of frequencies (1 to 80Hz, see [START_REF] Herrmann | Human EEG responses to 1-100 Hz flicker: Resonance phenomena in visual cortex and their potential correlation to cognitive phenomena[END_REF]) which in the context of Brain-Computer Interfaces (BCI) offers the possibility to present an extended number of classes with discriminative signal features. Hence, SSVEP-BCI are characterized by their high Information Transfer Rate (ITR) and classification accuracy. For instance, SSVEP-BCI based systems using classification algorithms based on Task-Related Component Analysis (TRCA) reported very high ITR (over 200 bits/minute) for a 40-target application (i.e., keyboard) thus offering a wide range of interactive commands [START_REF] Nakanishi | Enhancing detection of SSVEPs for a high-speed brain speller using task-related component analysis[END_REF]. Furthermore, the minimal training (e.g., using calibration-free classification methods based on Canonical Correlational Analysis (CCA), [START_REF] Lin | Frequency recognition based on canonical correlation analysis for SSVEP-Based BCIs[END_REF]) required by some of the SSVEP-BCI methods further established SSVEP paradigms as a reliable and efficient approach for BCI applications.

Nonetheless, current SSVEP-based BCIs have several pitfalls. Like many other paradigms used in BCI, SSVEP requires the user to direct his attention and gaze toward visual stimuli related to the command that corresponds to his intention. This constraint usually results in back-and-forth eye movements, and therefore divided attention, between the interface where the stimuli are displayed and the actual object affected by the commands. Another important concern regarding SSVEP-based BCI is the ocular fatigue induced by the prolonged presentation of flickering visual stimuli [START_REF] Cao | Objective evaluation of fatigue by EEG spectral analysis in steady-state visual evoked potential-based brain-computer interfaces[END_REF]. The strength (i.e., magnitude in terms of luminance/contrast modulations) of the visual flickers is positively correlated with the induced SSVEPs amplitude response [START_REF] Norcia | The steady state VEP in research[END_REF]. As a result, SSVEP paradigm designs have generally focused on maximizing SNR through the presentation of very salient, high-energy, visual stimuli. Furthermore, the amplitude of the SSVEP was found to be maximal when the stimuli flickers at a relatively low rate (ranging from 8 to 18Hz) [8] [9]. The repetitive presentation of visual information within such a low/intermediate frequency range is particularly straining on the eyes [START_REF] Cao | Objective evaluation of fatigue by EEG spectral analysis in steady-state visual evoked potential-based brain-computer interfaces[END_REF] and has been reportedly deemed annoying/frustrating by the users [START_REF] Ortner | An SSVEP BCI to control a hand orthosis for persons with tetraplegia[END_REF]. Another important consideration regarding the flickering frequency lies in the risk of epileptic seizures. Frequencies ranging from 4 to 20Hz present acute risks of triggering epileptic seizures [11] [12]. According to the prevalent theory of a limited pool of cognitive resources [START_REF] Kahneman | Attention and Effort[END_REF], sensory information whose lowlevel features are particularly salient drive bottom-up (i.e., object-oriented as the SSVEP stimuli) sensory and attentional processing at the expense of top-down (i.e., cognitive driven) processing strategies. In contexts where cognitive resources need to be allocated to attend to different aspects of a task, SSVEP stimuli may grasp attention away from the actual task at hand. Introducing such interference may be detrimental to the performance of the task. This is in clear opposition to the very aim of BCI applications directed at supporting individuals in the execution.

In order to make the SSVEP-based BCI accessible to a wide population and increase its comfort of usage, it is crucial to reduce the intrusiveness of visual stimulations. This study explores the reduction of stimulation amplitude depth for the design of SSVEP-based BCI applications as a potential solution to the aforementioned issues. To meet this goal, we compared classic flickers (100% magnitude depth) and low-amplitude (10% magnitude depth) flickers in terms of subjective rating, signal-to-noise ratio, and classification accuracy using two classification pipelines namely CCA, TRCA. Decreasing the stimulation amplitude depth reduces the contrast between the two alternating states of the stimuli (ranging from grey to white for the control high amplitude condition and ranging from grey to lighter grey for the low-amplitude condition) but also reduces their maximal luminance (high amplitude stimuli have a maximum luminance 10 times higher than low amplitude stimuli).

II. MATERIAL AND METHOD

A. Participants

Eight healthy volunteers (2 women, median age: 28 years, range: 24-29 years) with normal or correctedto-normal vision took part in this study. Participant 2 withdrew his participation reporting a high level of discomfort following the presentation of high amplitude stimuli (epilepsy screening questionnaire was negative). As a result, an additional participant (Subject 9) was recruited as a replacement to ensure a proper counterbalancing of experimental blocks presentation order. The participants did not report any of the exclusion criteria (neurological antecedents, being under psychoactive medication at the time of the study). The study was approved by the ethics committee of the University of Toulouse (CER approval number 2020-334) and was carried in accordance with the declaration of Helsinki. Participants gave informed written consent prior to the experiment.

B. Apparatus and stimuli

Visual rhythmic stimuli were presented on an LCD monitor (LG GN750 27-inch, 1ms IPS display, NVIDIA G-sync compatible), with a refresh rate of 240 frames per second, a resolution of 1024x768 pixel (width x height) and a luminance of 400 cd/m2. The visual rhythmic stimuli consisted in the presentation of series of squares (side of 200 pixels) at the centre of the screen whose colour state was alternating between grey and white at specific frequencies following a square function (see Figure 1). The choice of stimulus frequencies for the different classes was based upon 1) avoiding using a lower range of frequency to reduce eye strain [START_REF] Ortner | An SSVEP BCI to control a hand orthosis for persons with tetraplegia[END_REF] [6], 2) avoiding endogenous activity found in the alpha range (8-12Hz), 3) the restrictions relative to the refresh rate of the monitor (temporal resolution is limited by the integer number of frames that can be drawn every period). As a result, the following frequencies were selected: 13.33, 17.14, 21.81 and 24Hz. Each stimulus type (4 frequencies x 2 different magnitude depths) was presented 30 times for a total of 240 trials spread over four randomized blocks of 60 trials. Every trial lasted 6 seconds (4 seconds of stimulus presentation + 2 seconds of interstimulus interval, see Figure 1), therefore the experiment itself lasted for approximately 25 minutes including short breaks between blocks. The visual stimuli were generated using Psychopy functions [START_REF] Peirce | PsychoPy-Psychophysics software in Python[END_REF] running on Python 3.7. The visual flickers were presented against a grey background.

C. Experimental design and procedure

A four-class problem was used as a benchmark to compare the classification performance of low and highamplitude SSVEP stimuli. Participants were instructed to look at the centre of the screen where fixation crosses and stimuli were presented. Prior to the experiment, a series of visual stimuli flickering at the four different frequencies for both experimental conditions used in the experiment were presented to the participants to assess whether they experience eye strain or any discomfort related to the viewing of such stimuli by filling a self-reported questionnaire (stimuli were presented for 4 seconds after which participants filled-in the questionnaire, stimuli presentation order was randomized). The participants rated each stimulus on two Likert scales (ranging from 0 to 10) based on their subjective perception of how strongly the stimuli flickered and how comfortable the viewing of these stimuli was. The participants were then equipped with the EEG system. The capping procedure took from 20 to 30 minutes depending on the subject. As mentioned above, experimental blocks lasted for 25 minutes on average (including breaks between blocks) which amounts to a total participation time that did not exceed one hour.

D. EEG recording

EEG data were recorded from 32 active (Ag/AgCl) electrodes fitted in an elastic cap according to the 10-20 international system and connected to a LiveAmp amplifier (Brain Products, Munich, Germany). The ground electrode was placed at the Fpz electrode location and all electrodes were referenced to the FCz electrode. The electrode impedance was brought below 20 kΩ prior to recording. The data were acquired at a rate of 500Hz with an online digital band-pass filter ranging from 0.1 to 250Hz. At the onset of every stimulus presentation, an event trigger was generated by the stimulus presentation program and synchronized to the EEG data stream with Lab Streaming Layer (LSL) 1 .

E. EEG Processing and Analyses

The raw continuous EEG data have undergone an offline, high-pass filtering at 1Hz (zero phase, acausal, filter order: 1651, -6dB) and re-referenced to the average of all channels. Continuous EEG data were then epoched around event timestamps (-2 to 4 seconds epochs with flicker onset as time point 0). No channel was removed nor interpolated. Analyses were then carried out on the 4 seconds following stimulus onset. This study first characterized the magnitude of the SSVEP elicited for each frequency for both low and high amplitude conditions using Rhythmic Entrainment Source Separation (RESS, [START_REF] Cohen | Rhythmic entrainment source separation: Optimizing analyses of neural responses to rhythmic sensory stimulation[END_REF]). This comparison was performed on Signal-to-Noise Ratio (SNR) measures computed at RESS component level as recommended by [START_REF] Cohen | Rhythmic entrainment source separation: Optimizing analyses of neural responses to rhythmic sensory stimulation[END_REF]. RESS first computes channel-to-channel covariance matrices from narrow-band filtered data at stimulation frequencies (using a gaussian-shape filter of full width half maximum (FWHM) = .5Hz) and neighbouring frequencies (R matrices) (distance = 1Hz, neighbour filters FWHM = 1Hz). A generalized eigendecomposition is computed between the frequency stimulation and the average of neighbouring frequencies covariance matrices. The eigenvector with the largest eigenvalue is selected as the main RESS component. This latter component is then back-projected (acting as a spatial pattern) to the time series EEG data to maximize the SNR of SSVEP responses.

F. BCI Analyses

In order to assess the viability of using low amplitude flickering stimuli, the performance of two state-of-the-art classification pipelines were contrasted: one based on CCA and TRCA.

1) Canonical Correlation Analysis (CCA)

In recent years, CCA has become the gold standard method for the detection of frequency information in SSVEPbased BCI [START_REF] Lin | Frequency recognition based on canonical correlation analysis for SSVEP-Based BCIs[END_REF], [START_REF] Bin | An online multichannel SSVEP-based brain-computer interface using a canonical correlation analysis method[END_REF]. Indeed, CCA-based methods have demonstrated high classification performance (ITR of 267 bits/min for an online 40-target SSVEP speller [START_REF] Chen | High-speed spelling with a noninvasive brain-computer interface[END_REF]). One of the main advantages of the standard CCA is that it does not require any training or calibration data. Formally, it can be formulated as follows: let be X and Y two multivariate random variables. CCA looks for two weight vectors w x and w y that maximize the correlation between w T X X and w T y Y : arg max

(wx,wy) corr(w T x X, w T y Y)
This maximal correlation is the so-called canonical correlation.

In our context, X ∈ R Nc×Ns corresponds to a matrix of N c EEG channels by N s samples for a given epoch of data 1 https://github.com/sccn/labstreaminglayer and Y f ∈ R Nw×Ns corresponds to a matrix with the same length of X composed by N w artificial sine-cosine waves for the target frequency f, as well as their first harmonic. Then, one looks for the frequency f of the artificial sine-cosine waves Y f that maximizes the canonical correlation with our signal of interest X. We followed recommendations from [START_REF] Bin | An online multichannel SSVEP-based brain-computer interface using a canonical correlation analysis method[END_REF] regarding channel selection (O1, Oz, O2, P3, Pz, P4) and the number of harmonics (1) of the artificial sin-cosine waves to consider.

2) Task-Related Component Analysis (TRCA)

The second method, TRCA [START_REF] Tanaka | Task-related component analysis for functional neuroimaging and application to near-infrared spectroscopy data[END_REF] has been recently applied to SSVEP-BCI data recorded with both wet [START_REF] Nakanishi | Enhancing detection of SSVEPs for a high-speed brain speller using task-related component analysis[END_REF] and dry [START_REF] Xing | A high-speed SSVEP-based BCI using dry EEG electrodes[END_REF] EEG systems with remarkable results. The method has achieved around 90% of classification accuracy for 12 and 40 classes. While this approach requires the use of calibration data, it provides substantial improvements over other classification methods. For instance, [START_REF] Nakanishi | Enhancing detection of SSVEPs for a high-speed brain speller using task-related component analysis[END_REF] achieved an ITR value of 325.33 bits/min using a dry electrodes system. The TRCA method takes advantage of the filter bank analysis as described in [START_REF] Nakanishi | Enhancing detection of SSVEPs for a high-speed brain speller using task-related component analysis[END_REF], with the same parameters. Thus the following steps are repeated across the m bands obtained with the filter bank. In this study m = 5 sub-bands were set to m × 13 and 90Hz, following [START_REF] Nakanishi | Enhancing detection of SSVEPs for a high-speed brain speller using task-related component analysis[END_REF].

Filtering

First, the method uses calibration data to generate spatial pattern (a filter combining EEG channels) w (m) f for each class.

The spatial pattern for frequency f is obtained by solving the following optimization problem:

w (m) f = arg max w w T Sw w T Qw
Which can be solved by finding the leading eigenvector of the matrix Q -1 S, where Q is the covariance (∈ R Nc×Nc ) of the concatenated matrix of all trials with a stimulation frequency f and S is the sum of all covariances (∈ R Nc×Nc ) between separate trials for that same stimulation frequency. w T Sw is the inter-trials covariance of filtered data. w T Qw is the variance of the filtered data, it need to be bounded to obtained a finite solution for w (m) f .

Individual templates construction

The individual data templates

X (m) f 1, X (m) f 2, ..., X (m)
f N are obtained by averaging all the calibration data epochs (on the seconds following stimulus onset) for the corresponding frequency (i.e., class) and band of the filter bank. These calibration data are then discarded and not used for performance computation. After the calibration, any new data observation is divided into m sub-bands, and then the correlation coefficient between w T (m) X (m) and w T (m) X (m) is computed. The correlation scores for each sub-band are combined using a weighted sum (see [START_REF] Xing | A high-speed SSVEP-based BCI using dry EEG electrodes[END_REF] for more details). After repeating the process for all N frequencies, the maximal value across the N correlation coefficients is interpreted as the estimated frequency of the trial. The subjective rating of stimuli comfort was mainly affected by both amplitude reduction [F(1,7) = 93.444, p < .001, η 2 = 0.752] and the stimulus frequency [F(1,7) = 9.035, p < .001, η 2 = 0.066]. Moreover an interaction between the amplitude and the frequency of the stimulus was found to affect how comfortable the stimuli were rated by the participants [F(3,21) = 3.789, p = 0.026, η 2 = 0.026]. Posthoc t-tests (with a Holm correction applied) revealed that the low amplitude stimuli were rated as more comfortable (mean difference = 3.625) than the high amplitude stimuli [p < .001, d = 3.418]. Moreover, stimuli at higher frequencies were rated as more comfortable than stimulus flickering at lower frequencies (see Figure 2).

B. SSVEP SNR

The SSVEP SNR of the main RESS component for each stimulus frequency across both amplitude modulation conditions (low and high) are presented in Table I. A 2 × 4 repeated measures ANOVA was performed to investigate the effects of stimulus amplitude (low and high) and frequency A main effect of the amplitude reduction [F(1,7) = 17.421, p = 0.004, η 2 = .587] was observed on SSVEP SNR. In contrast, neither the stimulus frequency [F(3,21) = 0.611, p = 0.615, η 2 = .006] nor the interaction between stimulus amplitude reduction and frequency [F(1.4,9.6) = 0.383, p = 0.616, η 2 = 0.005] did not yield a significant effect on SSVEP SNR. Greenhouse-Geisser correction was applied when the sphericity assumption was violated. Post-hoc paired sample t-tests (with a Holm correction applied) revealed that low amplitude stimuli had a significantly lower SNR (mean difference = 20.5) than the high amplitude stimuli [p = .004, d = 1.476], see Figure 3.

C. BCI results

The results of the offline analyses performed on the whole stimulus presentation period are listed in Table II. For the high amplitude condition, the mean accuracy across all subjects reached 98.2% ± 2.1 with the CCA method and 99.8% ± 0.6 with TRCA. For the low amplitude condition, a mean accuracy of 62.8% ± 20.7 and 80.2% ± 15.8 were performed by the CCA and TRCA methods respectively. The minimal accuracy was observed in the low amplitude condition for subjects 1,8 and 9.

Figure 4 shows the accuracy of CCA and TRCA methods over different data lengths (0.5, 1, 2, 3 4 seconds). It can be observed that the TRCA method not only achieves higher classification accuracy for both the high and low amplitude conditions but also that the variance across subjects is also lower than CCA. Importantly, TRCA requires shorter data lengths than CCA to reach high classification accuracy and CCA performance continues to improve with longer data lengths whereas TRCA shows diminishing returns after 2 and 3 seconds of data for the high and low amplitude conditions. This result is of particular interest for the design of responsive BCI that operates on short time series. Overall, the TRCA method showed higher classification performance than the standard CCA approach for both high [t(7) = 2.531, p = 0.039, d = .895] and low amplitude [t(7) = 3.151, p = 0.016, d = 1.114] stimuli.

IV. DISCUSSION

The main motivation of this study was to assess the design of low amplitude SSVEP-based stimuli to improve user comfort and as a first step toward the implementation of transparent BCI. To this end, we have designed an experiment to compare SSVEP stimuli, presented at four different frequencies, under two magnitude levels: high (100%) and low (10%). As expected, the subjective rating scores revealed that low amplitude stimuli were perceived as more comfortable than the high amplitude stimuli. It is worth noting that higher frequencies were also deemed as more comfortable than lower frequencies. These results further suggest that decreasing the intensity but also increasing Is the increase of comfort with these low amplitude stimuli enough to counterbalance the frustration caused by longer stimulation and lower performances ? Nevertheless, it is worth noting that high classification accuracy was observed for some participants (100% accuracy for P6 and P7). The wide range of classification performance for the low amplitude condition hints toward interindividual differences. It is however unclear whether these differences in classification accuracy are related to individual differences in terms of stimulus perception sensitivity or nonlinear interindividual differences in SSVEP responsivity to various degrees of stimulation intensity. Moving forward, a systematic investigation of individual perception thresholds and a gradation of SSVEP responses at different levels of amplitude modulation would shed light on this pending question. Moreover, there is a need for a definition of subjectspecific breaking points of such a low amplitude approach. This information could be used to inform the definition of optimal amplitude selection to strike the right balance between SSVEP signal strength and user comfort.

Eventually, our results demonstrate that the TRCA algorithm outperformed CCA for both low and high SSVEP. This result is in line with recent studies reporting better classification performance using TRCA [19] [4]. TRCA seems to be an effective approach to deal with the reduced signal strength of low amplitude RVS and should therefore be considered as a classification method moving forward. In contrast to TRCA, CCA does not require to be trained on prior data. Such a calibration-free approach is particularly desirable from a user-experience perspective as it reduces the set-up time of the BCI. Nevertheless, calibration duration can be optimized with more interactive and engaging protocols to reduce the perceived constraints. Further investigations into classification approaches applied to low amplitude stimuli will shed light on the range of solutions available.

V. CONCLUSION AND PERSPECTIVES

This study was dedicated to assessing the relevance of low amplitude flickers for BCI purposes. Our results, collected over 8 participants, revealed that low amplitude visual stimuli elicited typical SSVEP responses while improving participants' subjective visual comfort. The encouraging classification performance achieved using the TRCA classification method further demonstrate the validity of using low-intensity flickers within the frame of SSVEP-BCI applications. Future work could extend the proposed approach to find the optimal compromise between global user experience and SSVEP signal strength following a subject-dependent and iterative approach. Moreover, the interaction between amplitude and other features of the flickering stimuli such as colour, duty cycle, patterns blur, or waveform shape may be studied to further reduce the intrusiveness of the visual stimuli. Eventually, the use of high flickering frequency could also present an alternative approach to promoting transparent and seamless brainmachine interaction for neuroergonomics applications [START_REF] Dehais | A Neuroergonomics approach to mental workload, engagement and human performance[END_REF].

Fig. 1 :

 1 Fig.1: Schematic representation of a single-trial SSVEP stimulus. Each trial is preceded by a 2 seconds (or 480 frames) interstimulus interval period during which a fixation cross was presented at the center of the screen. SSVEP stimuli are then presented for 4 seconds (or a total of 1920 frames). The SSVEP-eliciting stimuli is in practice the variation of transparency amplitude for a square presented at the center of the screen. The transparency follows a square waveform with a 50-50 duty cycle. The magnitude of the low amplitude flickers (red line) is set to be 10 times lower than in the high amplitude condition.

Fig. 2 :

 2 Fig.2: Distribution of mean visual comfort score for both amplitude modulation conditions (low and high) and each stimulus frequency(13.33, 17.14, 21.81, 24Hz)

Fig. 3 :

 3 Fig. 3: Grand average SSVEP Signal-to-Noise Ratio (SNR). A) Power spectrum of the different stimulus frequencies (13.33Hz -blue, 17.14Hz -orange, 21.81Hz -yellow, 24Hz -purple) for both amplitude modulation conditions (highdotted line, low -solid line). B) Mean time-course of SSVEP SNR over the length of stimulus presentation period. C) Topographical distribution of the RESS main component.

Fig. 4 :

 4 Fig. 4: Mean classification accuracy (with std. err.) over different data length (0.5, 1, 2, 3 4 seconds) for high (red solid line) and low (blue dashed line) amplitude stimuli using TRCA (square markers) and CCA (circle markers) methods.

TABLE I :

 I Main RESS component Signal-to-Noise Ratio (SNR) over the whole 4 seconds epochs for each stimulation frequency (13.33, 17.14, 21.81, 24Hz) and amplitude modulation (high and low).

TABLE II :

 II Comparison of classification accuracy obtained using the Canonical Correlation Analysis (CCA) and Task-Related Component Analysis (TRCA) methods for both amplitude modulation conditions.

ACKNOWLEDGEMENT

This work was funded by AID (Neurosynchrone project) and ANITI (Chair for Neuroadaptive Technology).