Morphological diversity, evolution and biogeography of early Pleistocene rabbits (Genus Oryctolagus)

The early Pleistocene is the key period for understanding the evolutionary history and palaeobiogeography of rabbits (Genus Oryctolagus). In western Europe, many species were endemic, making them a reliable indicator of the evolution of the terrestrial ecosystems in which many species have evolved. However, the morphological variability of rabbit species is still poorly understood and their phylogeny remains a subject of debate. Through both qualitative (morphological description) and quantitative (linear measurements and two‐dimensional geometric morphometrics) approaches, we address here the morphometric diversity of the third lower premolar (p3), a tooth commonly used to distinguish leporid species, in order to assess intra‐ and inter‐regional morphological variations in several early Pleistocene rabbit populations. Our results suggest that the different approaches are complementary and allow, on different levels, a full characterization of the p3 variability of early Pleistocene rabbits and to imply relations between populations. The size and shape variations of this tooth reflect the taxonomic and phylogenetic signals of the different species but were probably also significantly impacted by geographical position and local climatic conditions. In view of the great morphometric variability highlighted in this work, we suggest a careful reconsideration of certain dental criteria previously considered ‘diagnostic’ in the characterization of these species. However, the overall results allowed us to discuss the phylogeny of the genus Oryctolagus and to hypothesize the ecological requirements and different phases of the dispersal of taxa in western Europe, probably associated with global climate changes.

T H E quantification of morphological diversity is essential for achieving a better understanding of the evolutionary processes of species in time and space (e.g. Foote 1997; Hopkins & Gerber 2017). Studies of fossil specimens generally attempt to apprehend the factors responsible for morphometric and phenotypic variations (e.g. environmental and/or evolutionary event), their hypothetical phylogenetic links and the ecological characteristics that would make them valid and distinct species (Mayr 1963). This would prevent the amplification of specific differentiations (i.e. when anatomical characters differ between individuals) and the over complication of taxonomy. In absolute terms, these studies should be considered in the same comparative framework. However, the characterization of morphological diversity in a taxonomic group is often based on the use of distinct methods by the different authors, such as the description of morphological characters, traditional (i.e. using linear measurements) or geometric (i.e. landmarks and/or semi-landmarks) morphometrics. This has often hampered attempts to explore the interrelationships of different taxa, as well as their respective morphometric variability.
Descriptive approaches to morphology tend to focus on cladistic data sets in broad taxonomic scope studies (Lloyd 2016; Gerber 2019). They also allow us to grasp unique morphological characters. However, these approaches are less effective in addressing more tenuous changes in size, shape and/or proportion among individuals of the same genus or species. Conversely, traditional (i.e. linear measurements) or geometric (i.e. landmarks and/or semilandmarks) morphometrics can quantitatively characterize this kind of variation between specimens in a single study (Bookstein 1991;Zelditch et al. 2012). One of the strengths of these quantitative methods, and in particular the geometric morphometrics (GMM) now widely used in palaeontology, is identifying the processes that are at the core of these variations. Thus, none of these methods can necessarily be considered more effective since they are used in different circumstances. However, when we are interested in a particular taxon, it is often difficult to combine all of these different data from the scientific literature to assess and refine general macroevolutionary patterns. To solve these various issues, recent studies have compared these different qualitative and quantitative approaches to highlight whether patterns of morphological variation are congruent or contrasting (Hetherington et al. 2015;Hopkins 2017;Mongiardino Koch et al. 2017;Ferr on et al. 2020;Schaeffer et al. 2020). To this end, we explore the evolution of morphological diversity in rabbits of the genus Oryctolagus, the most abundant leporids in the Quaternary palaeontological and archaeological records from southwestern Europe. Over the past two decades, numerous old collection revisions and new material descriptions have allowed the taxonomy of Pleistocene taxa and their biochronological and palaeogeographical frameworks to be significantly improved. Despite this, the taxonomy and phylogeny of Pleistocene species still remain in flux and the subject of debate (  Ultimately, it would appear that all forms of the early Pleistocene coincide with specific periods and/or geographical areas (Fig. 1 To date, no inter-regional and inter-specific study has been conducted to analyse the diagnostic value of F I G . 1 . Palaeobiogeographic and temporal distribution of the genus Oryctolagus data in the late Pliocene and early-middle Pleistocene, based on published data. Squares represent the sites analysed in this study. Sites: 1, Layna, 2, Gorafe 2, 3, C ordoba (L opez- Mart ınez 1989Mart ınez , 2008 morphological characters in early Pleistocene rabbit species. Thus, it is quite conceivable that the differences observed between the fossil populations may be the result of chronological and/or geographical variations linked to environmental conditions, omitting the notion of intraand inter-specific variability. In addition, early Pleistocene rabbit fossils remain relatively rare, often fragmented and with a wide spatiotemporal distribution. Traditionally for the genus Oryctolagus, apart from studies on the quantification of variations in body size (e.g. Donard 1982;Sharples et al. 1996;Callou 2003;Davis & Moreno-Garc ıa 2007;Pelletier et al. 2015a, b;Pelletier 2018Pelletier , 2019Davis 2019), the vast majority of palaeontological works is based on the study of shape variations of the third lower premolar (p3). This tooth is widely used to distinguish leporid species and to reveal their phylogenetic relationships (Petter 1959(Petter , 1961Hibbard 1963;L opez-Mart ınez et al. 1976;Palacios 1976Palacios , 1996Palacios & L opez-Mart ınez 1980;L opez-Mart ınez 1980L opez-Mart ınez , 1989Donard 1982;Suchentrunk et al. 1994;Nocchi & Sala 1997a;Suchentrunk et al. 2000;De Marf a & Mein 2007;De Marf a 2008, 2009Angelone & Rook 2012;Angelone 2013;Pelletier et al. 2015a, b;Pelletier 2018Pelletier , 2019. However, the relevance of these dental criteria could be revised in relation to the data acquired on modern individuals. Indeed, it has been shown that some of these variations could be caused by selective processes such as adaptation to environmental conditions or regional gene flows (Suchentrunk et al. 2000). Moreover, it has also been demonstrated that latitude, temperature and precipitation (among other environmental parameters) are the main factors that drive the shape variation of p3 in the various modern populations of O. cuniculus (Pelletier 2019). This suggests the existence of a significant polymorphism and would call into question certain criteria that are generally used to identify rabbit species. In addition, the specific characterizations of fossil individuals using a description of these morphological characters could vary across authors, be subject to biases resulting from intrapopulation variations and potentially result in lack of consensus regarding the subjectivity of the criteria used to define a species.
Since it has been demonstrated that genetic divergence significantly impacts tooth morphology (e.g. Polly 2001;Caumul & Polly 2005;Ledevin et al. 2016;Cucchi et al. 2017Cucchi et al. , 2019, a more detailed analysis of dental shape using GMM could allow a distinction to be made between species or an assessment of taxonomic and phylogenetic signals. These approaches, now increasingly used in leporids (Patnaik 2002;De Marf a 2009;Pelletier et al. 2015a, b;Sawaura et al. 2018;Pelletier 2018Pelletier , 2019, could then enable the quantification of the variability of rabbit populations from the early Pleistocene and the recognition of how these could be related instead to intrinsic or extrinsic factors. Thus, in this work, we use three distinct approaches to characterize the morphological variability of the p3: (1) the description of dental characters; (2) the analysis of linear measurements; and (3) a two-dimensional (2D) GMM analysis of seven early Pleistocene populations, representing the three species described for this period (i.e. O. lacosti, O. valdarnensis and O. giberti). Thus, geometric morphometric characterization of the p3 will provide a fair and effective comparison with the two other methods more commonly used in studies of fossil leporids. It will also provide an opportunity to assess and discuss the relative effects of different factors (e.g. size, phylogeny, geographical position and climate) on morphological differentiation, in order to interpret the temporal, phylogenetic and ecological context of the evolution of early rabbits.

Sample
The analysed fossil sample includes 76 individuals from seven early Pleistocene populations representing the three species described for this period: O. lacosti, O. valdarnensis and O. giberti (Table 1)

Description and tooth measurements
The description of dental characters and the measurements of the p3 are based on the work of Donard (1982) (Fig. 2). Metric data were collected using a Dino-Lite Digital Microscope (AM7013MT). Our sample is described, compared and discussed with previously published data, including the type-species and affiliate specimens described for the early Pleistocene.

Geometric morphometric analysis
Data were collected from 2D digital photographs of the occlusal surface of each p3 using a standardized protocol (Pelletier 2019): teeth were positioned in a direction perpendicular to the focal plane, with parallax controlled for by symmetrically adjusting the lingual anteroconid and entoconid to limit biases resulting from the vestibulolingual obliquity of the tooth. Pictures were taken using the same Dino-Lite Digital Microscope. This 2D protocol involves six landmarks positioned on the inner edge of the enamel and corresponds to the points of maximum curvature (Fig. 2). Between landmarks 1 and 2, and 6 and 1, two curves, each composing 30 equidistant semilandmarks, were used to characterize the contours of the anteroflexid and the vestibular (or labial) and lingual anteroconids. As it has recently been shown that the morphological variations on the mesial portion of the tooth are directly impacted by different geographical and climatic factors, this protocol is relevant to identify what comes under taxonomic or extrinsic factors (Pelletier 2019). The hypoflexid was excluded from the analysis because although many authors report it as a trait of great phylogenetic importance, it does not appear to be related to taxonomy. For example, it has been shown that the number of crenulations increases with age and that the shape of these crenulations is correlated with their number, which is very variable within and among individuals (Bertonnier-Brouty 2019). Details of each specimen including catalogue numbers and repositories are provided in Appendix S1.
The landmark and semi-landmark coordinates were derived from digital photographs using tpsDig2 v.2.16 (Rohlf 2010a). All specimen coordinates were aligned using the generalized Procrustes superimposition to remove the position, scale and orientation parameters from the initial point configurations (Bookstein 1991). Unlike landmarks, semi-landmarks do not have an exact correspondence on the curvature of the enamel, and instead 'slide' along the line between adjacent points to minimize the sum of the Procrustes distances between each individual and the average shape (Perez et al. 2006;Gunz & Mitteroecker 2013). The generalized Procrustes superimposition of landmark and semi-landmark coordinates was conducted using tpsRelw v.1.49 (Rohlf 2010b).
To assess the order of magnitude of the morphological variation, as well as to identify whether the main source of this variation corresponds to chronology or taxonomy, we first compared the p3 morphology of early Pleistocene rabbits with that of modern populations (O. cuniculus). To do this, the shape differences between groups were estimated using a multivariate analysis of variance (MAN-OVA), with significant interaction (a = 5%) assumed to reflect the population and/or taxonomic differences. We then assessed the specific and/or population assignment accuracy by calculating the cross-validated correct classification percentages on the fossil sample, using a canonical analysis of variance (CVA). So as not affect the crossvalidation results, we reduced the dimensionality of our data set by keeping the values of the main components expressing 95% of the total variance before each canonical analysis (Kovarovic et al. 2011). However, given the limited number of p3 specimens in some deposits, we were forced to only keep the groups that had a sufficiently large number of individuals to perform this analysis (n > 5). We then focused the study on fossil populations in order to morphologically compare only early Pleistocene taxa with each other. Size differences were evaluated from log-transformed centroid sizes, using Kruskal-Wallis tests with an error threshold set at a = 5%. Pairwise comparisons of the populations were performed using multiple Wilcoxon rank tests according to these different categories. To control the false discovery rate, a multi-comparison correction was applied to the pvalues using the 'Benjamini-Hochberg' method (Benjamini & Hochberg 1995). The shape differences between these different groups were also estimated using a MAN-OVA, and the shape variation was visualized using principal component analysis (PCA) based on Procrustes coordinates. The phenotypic similarities between groups were calculated from Mahalanobis distances derived from canonical variates and visualized with a neighbour-joining network (Mahalanobis 1936). Finally, allometry was assessed using multivariate regressions of shape variables on the log-transformed centroid sizes. All morphometric statistics were performed with Rstudio v.
90°with some small undulations also present in specimen FSL-211901 ( Fig. 3A). However, the anteroflexid is deeper for both Saint-Vallier individuals than for Perrier one. Finally, the paraflexid is well marked, with cement visible on the occlusal surface.
Oryctolagus valdarnensis: The population of Pirro Nord 13 (n = 43) includes p3 specimens with relatively variable morphological characters (Fig. 3B). With the exception of one individual whose vestibular anteroconid is larger than the lingual one, and five individuals for whom both anteroconids have a relatively similar size and shape, most p3 specimens have a larger lingual anteroconid than the vestibular one (n = 37). The vestibular anteroconid is generally ellipsoidal for most the individuals and longer than the lingual one (n = 30) but may also be rather rounded (n = 13). However, the lingual anteroconid is generally globular (n = 19) or even rounded (n = 13), though rarely subsquare (n = 6) or subtriangular (n = 5). The depth of the anteroflexid greatly varies in this population since it can be from very deep (n = 25; e.g. Fig. 3B, UNIFE P13-USC3-QB2_44) to moderately deep (n = 12), as observed in the O. lacosti specimens from Saint-Vallier or generally in the extant rabbit O. cuniculus (Fig. 3D). Some specimens show a relatively short anteroflexid (n = 6; e.g. Fig. 3B, UNIFE P13-USD16-QB3_7) as observed in O. lacosti from Perrier or in O. giberti (Fig. 3C). The depth of the paraflexid, although always present, also varies significantly among individuals. It is often deep (n = 25) with cement observable on the occlusal surface in most cases (e.g. Fig. 3B, UNIFE P13-USD16-QB3_7, P13-USD15-QB4_15), but it can also only be marked by a slight sinuosity on the lingual edge of the tooth (n = 12; e.g. Fig. 3B, UNIFE P13-USC3-QB2_44, P13-USC8-QB3_11), or even be virtually absent (n = 6). The protoflexid is generally V-shaped with an angle of 90°or less (n = 38). Only five individuals present an obtusely angled protoflexid. Finally, the protoflexid has small undulations on the vestibular edge in one third of cases (n = 14), as observed in O. lacosti (Fig. 3A).
Oryctolagus giberti: The p3 of X abia II (n = 2) have an ellipsoidal vestibular anteroconid and a wider and subtriangular lingual anteroconid (Fig. 3C). The anteroflexid is deeper for MUPREVA 115.165, 769 than for 115.165, 772, but both individuals have a paraflexid that has a slight sinuosity. The protoflexid marks an angle of around 90°, although it is more acute for 115.165, 772 and rather obtuse for 115.165, 769. At the Cueva Victoria, p3 instead have an ellipsoidal vestibular anteroconid and a wider rounded lingual anteroconid. The anteroflexid is generally quite shallow and the paraflexid is marked by a slight sinuosity. The protoflexid is rather smooth with an angle of 90°or less. Only one individual (e.g. Fig. 3C, a previously undescribed p3) is characterized by its globular anteroconids of similar sizes, a clearly deeper anteroflexid and a protoflexid with small undulations. In the Boisde-Riquet deposit, p3 have either an ellipsoidal (n = 9) or rounded (n = 6) vestibular anteroconid. The lingual anteroconid is rather rounded (n = 13), but two specimens have a subtriangular one. Generally, the lingual anteroconid is either larger (n = 8) or the same size as the vestibular one. The anteroconids are separated by an anteroflexid of very variable depth: four showed a very shallow depth anteroflexid, six were a little deeper and four could be considered deep. However, the paraflexid is always slightly marked. The protoflexid is generally U-shaped with an angle of 90°or less (n = 11) and only four have an obtuse angle. At the Vallonnet cave, the vestibular anteroconid is predominantly ellipsoidal (n = 5) although it is also rather rounded in three individuals. The lingual anteroconid is quite variable: when it is rounded (n = 5) or ellipsoidal (n = 1), its size is comparable to or wider than the vestibular anteroconid, while those which have a subtriangular shape (n = 2) are smaller than the vestibular anteroconid. Concerning the anteroflexid, its depth varies between shallow (n = 5) and very shallow (n = 3). The paraflexid is very slightly marked in most specimens (n = 6). Finally, the protoflexid almost always has an angle of 90°or more; a single specimen was found with an angle is less than 90°.
Oryctolagus cuniculus: The individuals represented here were collected from the Balaruc I and Escale cave deposits (Fig. 3D). They are among the first representatives of O. cuniculus from western Europe. Generally, these rabbits have significantly deeper anteroflexids than what we have previously seen in early Pleistocene individuals. Both anteroconids are large, elongated mesiodistally and quite symmetrical. The depths of protoflexids and hypoflexids, as well as the number of folds of the enamel along the latter, are extremely variable criteria and do not allow trends to be identified within the populations analysed.

Conventional tooth measurements
In terms of size, we note a clear difference between O. lacosti from Perrier (n = 1) and Saint-Vallier (n = 2) in central-eastern France, with a larger p3, and p3 from Montouss e 5 (n = 2) in south-western France (Fig. 4). A significant size difference is also present in Italy between O. valdarnensis from Pirro Nord 13 (n = 43) in the province of Foggia, with a smaller p3, and those from Upper Valdarno (n = 2) in Tuscany. Although there are some overlaps between the rabbits of Pirro Nord 13, Montouss e 5, Vallonnet, Bois-de-Riquet, Cueva Victoria and X abia II, it seems that in general O. lacosti and O. valdarnensis are the taxa with the largest p3. Conversely, the populations of south-eastern Spain affiliated to O. giberti (Barranco Le on 5, Cueva Victoria, Fuente Nueva 3, Hu escar I, Quibas, X abia II) have the smallest teeth. In general, these rabbits have a smaller p3 than those of the same species in south-eastern France (Vallonnet cave, Bois-de-Riquet). In addition, the members of the population from Sima del Elefante (Castile and Le on, north-western Spain) have a larger p3. Finally, we note that whatever species considered, there is a relatively large range of values in the dimensions of the p3, with significant overlaps. However, although some samples are relatively small, it seems that there are significant size differences in the same species depending on their geographical location.

Two-dimensional tooth morphology
For each p3, the factorial MANOVA found significant differences in shape among the different groups (all p << 0.05; Table 2). The discriminant model found correct chronological (i.e. early Pleistocene vs modern populations) and specific identification for 91% and 95% of the classifications after cross-validation, respectively (Table 2). This shows that this model can accurately distinguish an isolated p3 from the archaeological record in an early Pleistocene context. However for each site, the discriminating model only classified 69% of the p3 specimens. On the other hand, a clear taxonomic signal segregated the populations into three groups (Fig. 5): (1) the four modern populations belonging to O. cuniculus (i.e. Las Lomas, Navarre, Santar em and Tour du Valat); (2) the Boisde-Riquet, Cueva Victoria and Vallonnet populations related to O. giberti; and (3) the population of Pirro Nord 13 affiliated to O. valdarnensis. Although the differences are tenuous, the visualization of the shape changes along canonical axes 1 and 2 showed that in modern populations, the anteroflexid is deeper, vestibular and lingual anteroconids are symmetrical and the paraflexid is very weak compared to the early Pleistocene species. For the p3 specimens from Pirro Nord 13, the anteroflexid is shallower and the anteroconids are asymmetrical: the vestibular anteroconid is rectangular and more elongated mesially, while the lingual anteroconid is rounded. This asymmetry is explained by the presence of a very marked  paraflexid. Finally, the p3 of O. giberti from Bois-de-Riquet, Cueva Victoria and Vallonnet are closer to those from Pirro Nord 13, particularly concerning the shape of the mesial part of the tooth (CV1). However, the paraflexid remains present but less marked (CV2).
In the fossil populations investigated, a Kruskal-Wallis test produced significant differences in the centroid size of the p3 between the different rabbit species (p = 1.72 9 10 À6 ; Fig. 6A) as well as the different fossil populations (p = 5.14 9 10 À5 ; Fig. 6B). Although the Perrier, Saint-Vallier and X abia II deposits only yielded a relatively small number of individuals (all n < 3), the centroid size of the p3 of O. lacosti was the largest, while the centroid size of X abia II rabbits (O. giberti) was the smallest (Fig. 6). The populations of O. giberti did not differ significantly from one another (Table 3) but their p3 were always smaller than the rabbits of Pirro Nord 13, assigned to O. valdarnensis. In addition, MANOVA (with p = 7.79 9 10 À4 ) showed a clear structuring of p3 shape variation (Fig. 7). Individuals from Perrier and Saint-Vallier had a similar shape and showed greater affinity with the population from Pirro Nord 13. The populations of south-eastern France (Bois-de-Riquet and Vallonnet) and Spain (Cueva Victoria and X abia) showed a close dental shape. Variations in p3 shape in early Pleistocene rabbits correspond to a deepening/shallowing of the anteroflexid (PC1) associated with the enlargement/ narrowing of the lingual and vestibular anteroconids (PC2; Fig. 8). The population of Pirro Nord 13 had a relatively deeper anteroflexid, a marked paraflexid, as well as a large lingual anteroconid, rounded and projecting towards the inside of the mandible. In the O. giberti populations, the anteroflexid and the paraflexid were more reduced and the anteroconids are proportionally smaller and less rounded. Finally, O. lacosti from Perrier and Saint-Vallier had a shape that was closer to the average shape of the sample analysed, with a paraflexid of intermediate size, as well as two large anteroconids of the same size. However, despite the large size differences previously observed in the centroid size of p3 (Fig. 6), the allometric component indicated no relationship between size and shape (p = 0.06).

Evolution and morphological variability of early Pleistocene rabbits
Overall, the morphometric descriptions and GMM analyses presented here have allowed, on different levels, a more optimal characterization of the morphological variability of fossil rabbit p3 in various early Pleistocene populations. Firstly, the p3 of O. lacosti has been described as presenting two large anteroconids with a more elongated and ellipsoidal vestibular anteroconid than the lingual anteroconid, which is wider. However, the depth of the anteroflexid seems to vary more than previously described by De Marf a & Mein (2007) (Fig. 3); it is relatively shallow in specimens from Perrier (FSL-211647; this work) and Montouss e 5 (F-120 and F171; Chaline et al. 2000, p. 101), while it is deeper in Saint-Vallier (FSL-211901 and FSL-495764, this work). In this species, the paraflexid is always present, although it is significantly more marked in Saint-Vallier than in the two other above-mentioned Gelasian deposits. Thus, the overall p3 conformation of Perrier is clearly closer to the individuals from Montouss e 5 than Saint-Vallier. In the Pirro Nord 13 deposit, despite significant inter-individual morphological variability, the p3 morphology coincides more with those of the few individuals of O. valdarnensis described in Valdarno and Pirro Nord 10 (Angelone & Rook 2012; Angelone 2013). The vestibular anteroconid is generally more elongated than the lingual anteroconid, which instead has a more globular shape. This is notably due to the presence of a very marked paraflexid in the majority of individuals, although it is shallower in a few individuals at Pirro Nord 10 (Angelone 2013 Indeed, although there may be some individual variations, the anteroflexid is generally shallow in all O. giberti populations. In this species, the p3 appears to have two anteroconids with similar size and shape, particularly in the French specimens from Bois-de-Riquet and Vallonnet, while the lingual anteroconid tends to be wider in the Spanish deposits. In Bois-de-Riquet and Vallonnet, the protoflexid is U-shaped and has an angle of 90°or less, while it is V-shaped at Cueva Victoria and X abia II. The paraflexid almost always has a slight sinuosity. Generally, the p3 also appears to be less elongated mesiodistally than in other species. These observations slightly differ from the description proposed by De Marf a (2008), which indicated in particular the presence of a deep anteroflexid in the holotype of O. giberti. Thus, the observations made in our sample are closer to the description proposed by L opez Jim enez et al. (2020) for the Cueva Negra rabbits (i.e. two large anteroconids with a similar shape and size, separated by a shallow anteroflexid and a well-marked paraflexid), with the exception of the protoflexid shape, which has obtuse angle. However, we agree with De Marf a (2008) that O. giberti can be distinguished from modern O. cuniculus, which generally has large symmetrical anteroconids, a deep anteroflexid and a more attenuated paraflexid. Thus, when we compare the different populations and species with each other, the descriptions of dental characters allow us to mainly identify the general tendencies within a population or group, rather than the characters that are strictly exclusive to each group. Indeed, although our observations are quite consistent at the species level with the descriptions of their respective lectotypes and holotypes (e.g. the depth of the anteroflexid and paraflexid or the general shape of the anteroconids), we also noted the presence of a relatively significant intra-and inter-population polymorphism which induces variations in the descriptions of the different authors. For example, De Marf a & Mein (2007) and Angelone & Rook (2012) described a p3 with a deep anteroflexid in O. giberti and O. valdarnensis respectively, but this depth is actually very variable, not only at the species level, but also within a given population. This could unfortunately result in the misidentification of an isolated tooth or a small sample. However, despite the reduced number of specimens for these periods, it is possible to identify trends in metric variations in p3 size related to the geographical position of the populations (Fig. 4). Indeed, O. lacosti has significantly larger dental dimensions in Perrier and Saint-Vallier, located in central-eastern France, than in Montouss e 5, located at a lower latitude in south-western France. The same observation was made in O. valdarnensis specimens, in which individuals from Valdarno in northern Italy also had a larger p3 than those individuals from Pirro Nord 13 in the south-eastern part of the country. Finally, the O. giberti from southern France are generally larger than those from south-eastern Spain. This can be explained by the fact that in western Europe, there is a strong correlation between the body size of rabbits (i.e. the size of dental and bone elements) and both geographical location and local environmental conditions France. Therefore, the size of dental and bone elements cannot be considered to be a reliable criterion for specific identification, since it can vary considerably and be directly related to local conditions. We also noted significant variability in the measurements taken by the different authors of the p3 of the same individuals, up to 0.3 mm of difference in the mesio-distal diameter (MDD) and vestibulo-lingual diameter (VLD). We noted, for example, a significant difference in the dental measure- Our 2D GMM study also provides additional information on morphological descriptions and variations in dental measurements, which revealed numerous overlaps between the different populations. First, it allows us to support the validity of the two species present during the Calabrian: O. giberti and O. valdarnensis. As had been previously hypothesized, our results confirm that all the rabbits that occupied the Iberian Peninsula (De Marf a 2008, 2009) (2012) who considered that all the early Pleistocene Oryctolagus remains collected thus far in Italy belong to O. valdarnensis. Morphologically, the analysis supports the hypothesis that the two Calabrian species are mainly distinguished by the depths of the anteroflexid and paraflexid, significantly deeper and more marked in O. valdarnensis than in O. giberti, leading to a notable dissymmetry of the two anteroconids. However, most of the variance is observed on axis 1 (26%) and is expressed on the intra-specific scale, thereby proving the great morphological variability in the different fossil populations in time and space (Fig. 8). This particularly suggests that there is a greater variation in the depth of the anteroflexid in each species, as well as a mesiodistal narrowing or enlargement of the occlusal surface, effectively impacting the size and shape of the anteroconids. As such, the recent study on different modern O. cuniculus populations from western Europe showed, for example, that the anteroflexid was shorter in south-western populations on the Iberian Peninsula (evolving in an environment in which the average annual temperatures are higher, average annual precipitation is lower and with a Mediterranean scrub-type vegetation) than in the populations of northern Spain or southern France (evolving in a more oceanic climate and with more precipitation and lower temperature, under a 'temperate forest' type vegetation) (Pelletier 2019). Thus, when we take into account the notion of intra-specific or intrapopulation variability, as well as the likely global/local environmental impact on this variation, our results show a degree of subjectivity in some of the qualitative morphological criteria for the description or redefinition of a species, particularly on the mesial part of the p3 (i.e. anteroflexid, paraflexid, anteroconids).
As expected, our discriminating analysis reveals that early Pleistocene species can be considered to be distinct from the extant species O. cuniculus. It is very important to take into account that this model can accurately distinguish an isolated p3 from the archaeological record in an early Pleistocene context, given both the great morphological variability and eco-ethology of the modern rabbit. Indeed, this species is a burrowing animal capable of producing complex and substantial underground warrens (Biadi & Le Gall 1993) and entering an archaeological or a palaeontological site naturally and/or accidentally (Pelletier et al. 2016(Pelletier et al. , 2020. This particularity increases the possibility of discovering its remains in the fossil record and poses questions about the contemporaneity of these remains with other assemblage components (Pelletier et al. 2015b(Pelletier et al. , 2016(Pelletier et al. , 2017(Pelletier et al. , 2020. Thus, in the case of old excavations and a lack of taphonomic context, our results show that the dental morphology of extinct species can be useful and complementary to taphonomic analyses in order to ensure the integrity of the archaeostratigraphy and the reliability of the studied material, as well as for biostratigraphic purposes. Nevertheless, the major concern of our study is to highlight the validity of the O. lacosti species. In this respect, Angelone & Rook (2012)  . Thus, only a few isolated remains have been distributed between the Perrier, Saint-Vallier or Montouss e 5 deposits, making it difficult to establish distinctive criteria for this species. In light of our results and both the qualitative and quantitative data from the literature, we are currently unable to make a clear distinction between O. valdarnensis and O. lacosti. However, we have noted that a few isolated individuals such as the individuals in Perrier and Valdarno were quite similar, whereas the morphometric characteristics of the Perrier and Saint-Vallier rabbits highlighted in the GMM study are relatively close to those observed at Pirro Nord 13 (Fig. 7). In any event, the individuals from Perrier and Saint-Vallier share significantly and statistically more morphological similarities with O. valdarnensis than with other species of rabbits. Consequently, it is entirely reasonable to envisage: (1) a phylogenetic succession between O. lacosti and O. valdarnensis, in which the latter would then have become the endemic form on the Italian Peninsula in the second half of the early Pleistocene, as already hypothesized (Angelone 2013); or (2) a synonymy between the two forms in which O. valdarnensis could be a more recent form than O. lacosti and these two forms should be considered to be geo-chronospecies, as previously proposed (Pelletier 2018). Moreover, due to the limited amount of data (individuals, deposits) widely dispersed across western Europe during this period, we do not have a sufficient number of elements to suggest any scenario. More O. lacosti material from Gelasian contexts would be needed to extend this debate. This would allow us to better evaluate the chronology of the Oryctolagus from this period to establish whether or not the French and Italian deposits were coeval. Indeed, the last O. lacosti from France are reported from around 2.1-2.0 Ma in Saint-Vallier while the first O. valdarnensis from Italy have been placed around 2.1 Ma in Montagnola Senese. However, these relative datings (i.e. not radiometric) are not accurate enough to establish whether or not they really overlap. However, despite the obvious phylogenetic proximity between the analysed individuals of O. lacosti and O. valdarnensis, we propose at this stage that these two taxa should be retained, representing only geo-chronologically these two rabbit forms in southern France and Italy, respectively.

Palaeoecology and palaeobiogeography of rabbits during the early Pleistocene
As well as understanding more about the phylogenetic relationships between the different rabbit species from the early Pleistocene, our study also makes it possible to discuss their diffusion dynamics in western Europe. Indeed, the expansion or isolation events of rabbit populations could have been conditioned by various ecological factors. Although the palaeoecology of the different early Pleistocene rabbit species has not been clearly established, their ranges were very close to that of the Pleistocene O. cuniculus (Pelletier 2018). Thus they could have been conditioned by the same biotic and abiotic factors. This would tend to assume that the Oryctolagus species shared, at least partially, similar ecological requirements to the extant European rabbit. These requirements could also explain why the genus Oryctolagus naturally spread only between the Iberian Peninsula, southern France and the Italian Peninsula during the Pleistocene, and not to northern France, central Europe or the Balkan Peninsula. Environmental components such as aridity, altitude and soil conditions are fundamental for the development of the extant European rabbit, mainly for the establishment of their burrows and for the maintenance of populations. Consequently, the European rabbit prefers loose and welldrained soils on flat or sloped terrain, such as sand dunes or grasslands, with low vegetation comprising bushes and shrubs at altitudes of up to 1000 m (Biadi & Le Gall 1993). It therefore does not colonize arid, overly dry or frozen (i.e. permafrost), mountainous, wetland or large forest areas. Thus, both topography and fluvial geodynamics are factors that have clearly conditioned the population migrations in western Europe (Fig. 9). In the south of the Iberian Peninsula, the three main rivers, the Tagus, the Guadiana and the Guadalquivir, have probably conditioned rabbit populations along the Mediterranean coast. Moreover, the mountains of the Iberian System in the centre of the Peninsula and the Sierra Nevada in the south, could have been an obstacle to an east-west circulation and, towards the north, the Ebro and the Pyrenees mountains would also have limited the dispersal of these species. Thus, the palaeontological data available for this period are in agreement with the genetic data, which have identified a refuge area in southern Spain and along the Mediterranean coast (Branco et al. 2002). The current lack of rabbit remains from the early Pleistocene in south-western France allows us to consider a single route of migration to the west of the Pyrenees, via Catalonia, between France and Iberia. The Garonne, the Tarn, the Loire and the Massif Central then preferentially orientated migrations towards the south-east. In south-eastern France, rabbits appear to have prospered on a narrow region between the Rhône delta, the Rhône Valley and the Southern Alps. Thus, the narrow passage between the Alpine chain and the Mediterranean Sea was the only possible route taken by Oryctolagus to reach Italy. Rabbits were then able to extend their distribution towards the south of the Italian Peninsula. The lack of palaeontological data thus far beyond the Northern Apennines chain suggests that these mountains have limited the diffusion of populations to the north and instead favoured a diffusion corridor along the west coast.
Although these geographical constraints may initially have guided the preferential migration routes of early Pleistocene rabbits, it is also likely that local and/or global environmental variations, such as temperature or precipitation, may have affected these dispersions. It has already been assumed that the climatic oscillations of the middle and late Pleistocene considerably influenced the geographical distribution of the extant rabbit in western Europe (Pelletier 2018). Indeed, although O. cuniculus is a species that is more adapted to a Mediterranean-type climate, it has adapted relatively well to brief or slight Pleistocene climatic variations. However, during sudden or extremely unfavourable climatic events, the rabbit occasionally disappeared from certain regions during the middle and late Pleistocene. Generally, it is during these critical climatic events that the fragmentation or isolation of continental biotopes can occur (Blondel 1986). These environment segmentations have regularly been observed during the Pleistocene glacial periods, generating veritable 'refuge areas' for different animal species (e.g. Hewitt 1996Hewitt , 2004Sommer & Nadachowski 2006;Valensi 2009). The isolation of populations in these refuge areas often leads to founder events responsible for speciation processes. It has also been shown that the rates of speciation and extinction significantly increase during these critical phases, highlighting a close link between ecological and evolutionary changes (G omez Cano et al. 2013). During climatic favourable phases, these processes can be followed by expansion(s) and/or recolonization(s) that have a significant impact at the genetic and phenotypic level (e.g. Mayr 1942;Hewitt 2000).
According to Nocchi & Sala (1997b, p. 183), all Oryctolagus had a common western European origin, which can be recognized either in O. laynensis or in an unknown progenitor form. As previously mentioned (see Fig. 1), O. laynensis seems to have been confined to the Iberian Peninsula until the end of the Piacenzian or perhaps even until the beginning of the Gelasian, and was known to be generally associated with an arid, warm subtropical savannah-type fauna (L opez-Mart ınez 2008). In southern France, O. lacosti were only present in the region from the Gelasian and evolved under a relatively humid but not particularly cold climate in a mosaic landscape comprising a mixture of steppe and open wooded areas (Gu erin et al. 2004). Thus, an allopatric speciation seems possible between these two Oryctolagus lineages, which are quite distinct ecologically and morphologically (L opez-Mart ınez 2008) and could have occurred following their geographical segregation by the Pyrenean chain.
These speciations, not necessarily synchronous, could be consistent with the sudden increases in aridity during the Pliocene (around 4.0, 3.6, 3.2 or 2.8 Ma), which led to significant changes in the small mammal communities of western Europe (Agust ı et al. 2001). The Gelasian is then characterized by a long period of climatic stability, favouring the dispersion of many mammals in southern Europe (Azzaroli et al. 1988;O'Regan et al. 2011;Bellucci et al. 2014). Although the lack of data in the Iberian Peninsula does not allow us to understand the evolution of O. laynensis during this period (i.e. maintenance, new speciation or disappearance), it was probably favourable for the maintenance and dispersal of O. lacosti. It was also during this period that the first Oryctolagus reached Italy (around 2.1 Ma). Whether O. lacosti and O. valdarnensis are synonymous or two distinct species, the endemic nature of O. valdarnensis, due to an isolation phase, could have resulted from the cold climatic event clearly marked at marine isotope stage (MIS) 82 (Popescu et al. 2010). It is also probable that after this climatic event Oryctolagus disappeared from southern France and was isolated in the peninsulas.
The Calabrian subsequently shows several major climatic changes (Agust ı et al. 2010). Thus, it is probably F I G . 9 . Main topographic and fluvial constraints that conditioned the dispersal routes of the early Pleistocene Oryctolagus in western Europe. Arrows represent the preferred migration routes; dots represent the same sites and species as in Figure 1. one of these major changes, marked by a general drop in temperature and an increase in aridity around 1.9-1.7 Ma (Agust ı et al. 2009;Bertini 2013), which, in turn, impacted the evolutionary and biogeographic dynamics of rabbits, isolated in the Iberian and Italian peninsulas during this period. In fact, it seems that it was after this climatic event that the first O. giberti appeared in Spain. This species was abundant in the south-east of the Iberian Peninsula, where the X abia II population could represent the oldest evidence of this species in Europe. The palaeontological occurrences would indicate a diffusion of the species which would mainly have occurred towards the north, along the Mediterranean coast. Around 1.4-1.2 Ma, this phase is characterized by a sharp increase in temperature and precipitation (Agust ı et al. 2009) followed by relative climatic stability with warm conditions (between 1.2 and 0.9 Ma; Altolaguirre et al. 2019). This could explain the rapid expansion of rabbits towards the north of the Iberian Peninsula and the south-east of France. In addition, we found that the O. giberti populations of south-eastern Spain had a smaller p3 size than those of south-eastern France (Figs 4 and 6), which could indicate that the size of rabbits during this period was influenced by the same latitudinal variations as modern populations (see above). However, locations with relatively similar latitudes can also have quite different climates (Peel et al. 2007), and thus present contrasting results. In modern populations, a peculiarity has been reported concerning the populations of Navarre (northern Spain), latitudinally lower than the populations of southern France but with equivalent or even slightly higher average size values (Callou 2003;Pelletier 2018Pelletier , 2019. In this case, the influence of the local environment overrides the effect of latitude. In fact, the populations of Navarre evolved under a climate subject to oceanic and mountainous influences, unlike the populations of southern France that live in a Mediterranean climate, where the annual average temperatures are relatively higher with lower pluviometry. Thus, climatic factors could have strongly influenced the population of Sima del Elefante in northern Spain, currently undergoing the same environmental conditions as in Navarre and with dental dimensions larger than those of the populations of southern France in the Vallonnet cave and Bois-de-Riquet, yet located at higher latitudes. It is also important to weigh these results against the micromammal associations of these deposits, which are not strictly coeval. For example, the micromammals reflect a generally temperate climate in Sima del Elefante (Rodr ıguez et al. 2011;Benn asar et al. 2016), while a cooler and drier climate is recorded in the Vallonnet cave (Montuire & Desclaux 1997; Moull e 2012), or a relatively more humid climate in Bois-de-Riquet (Lozano-Fern andez et al. 2019). Thus, it is crucial that geographical data are taken into account in the palaeontological studies of rabbits because this can provide additional elements of discussion that could enable a better reconstitution of the microclimate variations and local environments.
As we have noted, the period between 1.2 and 0.9 Ma was favourable to the expansion of the Spanish O. giberti populations towards south-eastern France. However, another major event that led to the disappearance of O. valdarnensis from the Italian fossil record around 1.2 Ma may correspond to the maximum cold peak of MIS 36 (Joannin 2007, Altolaguirre et al. 2019. This event was marked by a general change in the composition of forests and a major renewal of the mammal community in Europe (Head & Gibbard 2005; Magri & Palombo 2013). However, it would be important to include more material from O. valdarnensis in our study to gain a better insight into the impact of climatic conditions on migration and evolution dynamics across the entire Italian Peninsula, as is the case in Spain and France for O. giberti. Nevertheless, although body size gradually increases with latitude, a different impact has been noted between O. giberti and O. valdarnensis. Indeed, for similar latitudes, we have found that Italian rabbits are comparatively larger than individuals from south-eastern Spain and south-eastern France. This could be partially explained by the rather temperate and dry environments (Pavia et al. 2012;Combourieu-Nebout et al. 2015). However, the differential body size variations between the two peninsulas could also be an argument to suggest that these rabbits belong to two distinct phyletic groups. As such, despite the impact of environmental conditions on the body size of individuals, the allometric component was not significant in our GMM analysis. This suggests that the morphological differences between these two species are not only related to the variation in the body size of individuals, but could also reflect a taxonomic and phylogenetic signal. Finally, the transition between the early and middle Pleistocene is also marked by significant global cooling, causing deep mutations, particularly in the biomes of the northern hemisphere (Head & Gibbard, 2005;Joannin 2007). This cooling led to the fragmentation of landscapes and the isolation of various Villafranchian species in more favourable areas (Bennet et al. 1991;Bonifay & Brugal 1996). The different regions of Europe then underwent deep faunal renewal (Magri & Palombo 2013) which, in rabbits, corresponds to the extinction of O. giberti and the appearance of the extant rabbit O. cuniculus.

CONCLUSION
Most palaeontological studies have used p3 morphology to explore taxonomic and phylogenetic relationships in different Oryctolagus species. Unfortunately, most early Pleistocene species have been described or redefined from a small sample of individuals (i.e. n < 5), which is unlikely to accurately reflect the morphological variability of a particular population or species. Intra-and interpopulational studies appeared to be an important and unprecedented attempt to identify these variations, which are essential to more accurately specifying the phylogeny and biogeography of the genus. Overall, the different methods of analysis used in this study (i.e. morphological description, linear measurements and 2D GMM) have proven to be relatively complementary and enabled, on different levels, a characterization of the morphological diversity of early Pleistocene rabbits and to make connections among populations.
By means of a precise dental nomenclature widely used in palaeontology, morphological description seems to be adequate for describing an isolated individual or a small group of individuals, and to make small-scale comparisons. However, when the samples are larger, the descriptions are often limited to an average shape or the recurrence of main criteria. This method does not adequately capture the overall intra-and inter-populational or intra-specific variations, which, given the large morphometric variability of rabbits, could be problematic for identifying the species from an isolated p3. From one or two direct measurements of the tooth, linear measurement analyses enables the p3 size of an individual or group of individuals to be quickly captured and compared with other populations. However, as the body size of rabbits is directly related to their geographical position and local environmental conditions, this does not allow the identification of a given species. Nevertheless, to a certain extent, it could provide some clues about environmental variations in the immediate vicinity of a deposit. Lastly, GMM allows a comparison to be made of all individuals in a sample in the same morphological space and to directly visualize the results of statistical analysis in terms of conformation difference. This makes it possible to more easily recognize character variations that might be linked to intrinsic (i.e. phylogenetic) or extrinsic (i.e. ecological) processes. However, this method requires a relatively large data set and the establishment of a precise protocol, which takes longer to implement. Thus, GMM appears to be more adequate in a multi-regional study and over a long chronology, as was the case in this work. In this sense, neither of these approaches is necessarily more effective, and their use should depend on the research question under consideration.
Firstly, this work allowed us to confirm the specific status of two of the species previously defined in early Pleistocene, namely, O. giberti in the Iberian Peninsula and south-eastern France (probably throughout the Calabrian; i.e. 1.8-0.8 Ma) and O. valdarnensis in the Italian Peninsula (between 2.1 and 1.2 Ma). This also allowed us to phyletically relate the Gelasian species O. lacosti and O. valdarnensis. Due to the limited material available for O. lacosti, as well as an imprecise chronological position (i.e. relative, not radiometric), it is currently impossible for us to state whether these two species coexisted in western Europe or whether O. valdarnensis immediately succeeded O. lacosti in time and space. Further research, including more Gelasian populations, should allow us to refine our palaeobiogeographical knowledge and our hypotheses on the phyletic relationships that unite these two species. However, our study highlights the great morphological variability of rabbits from this period, at both the inter-specific and intra-populational level. Thus, although these species had previously been adequately recognized and defined from morphological descriptions and linear measurements, our results suggest that local environmental conditions are the main factors behind the variation in p3 shape. In fact, in this relatively important geographical and chronological context, this leads us to temper certain so-called 'diagnostic' criteria traditionally used to identify species, such as anteroflexid and paraflexid depths or even the shape of the anteroconids, because they might reflect ecological parameters rather than taxonomical variations. In this sense, these results strongly suggest that GMM analysis of dental traits are relevant markers that could be used to document the diversity and distribution of fossil rabbit species. Finally, geographical and chronological distributions and their associations with other mammals indicate that the ecological requirements of early Pleistocene species should partially correspond to those of the extant European rabbit. Thus, we can discuss population dynamics in Europe that were conditioned by topography, fluvial geodynamics and climatic and environmental variations. New palaeontological discoveries and the extension of studies to other dental and skeletal remains should ultimately allow even better clarification of the phylogeny, palaeoecology and palaeobiogeography of early Pleistocene rabbits in western Europe.
from Pirro Nord. I am very grateful to Marie Matu, the reviewers, editors and publications officer for their constructive comments that greatly improved the quality of the manuscript. Finally, this work was financed by an A*MIDEX grant (No. ANR-11-IDEX-0001-02) from the French Government programme 'Investissements d'avenir', and I would like to warmly thank the University of Oulu's archaeology lab for its support.

SUPPORTING INFORMATION
Additional Supporting Information can be found online (https:// doi.org/10.1111/pala.12575): Appendix S1. List of specimens studied including their catalogue number and repository, as well as the raw linear measurement of the p3 used in this study. Abbreviations: MDD, mesiodistal diameter; VLD, vestibulo-lingual diameter.