# The contribution of functional traits to the understanding of palaeoenvironmental changes

Léa Terray, Emmanuelle Stoetzel, Anthony Herrel, Raphael Cornette

► To cite this version:

Léa Terray, Emmanuelle Stoetzel, Anthony Herrel, Raphael Cornette. The contribution of functional traits to the understanding of palaeoenvironmental changes. Biological Journal of the Linnean Society, 2021, 133 (4), pp.1110-1125. 10.1093/biolinnean/blab057. hal-03233705

HAL Id: hal-03233705

https://hal.science/hal-03233705

Submitted on 13 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
The contribution of functional traits to the understanding of palaeo-environmental

changes

Léa Terray<fnref n="1" />*, Emmanuelle Stoetzel<fnref n="2" />, Anthony Herrel<fnref n="3" />, Raphaël Cornette<fnref n="1" />

<footnote n="1">¹ Institut de Systématique, Evolution, Biodiversité (ISYEB), Muséum National d'Histoire Naturelle, CNRS, Sorbonne Université, EPHE, Université des Antilles, CP50, 57 rue Cuvier, 75005 Paris, France.</footnote>

<footnote n="2">² HNHP, CNRS, Muséum National d'Histoire Naturelle, Sorbonne Université, UPD, Musée de l'Homme, Palais de Chaillot, 17 place du Trocadéro, 75016 Paris, France.</footnote>

<footnote n="3">³ Mécanismes Adaptatifs et Evolution (MECADEV), Muséum national d'Histoire naturelle, CNRS,
55 rue Buffon, Paris, France</footnote>

<footnote>* Corresponding author (mail: lea.terray@mnhn.fr)</footnote>
Abstract

Performance traits interact with the environment and are consequently relevant ecological
indicators. In particular, performance traits implicated in feeding may be powerful ecological
proxies are they are directly related to diet. However, performance traits have only rarely
been used to better understand palaeo-environmental variation. In the present study we
evaluate the usefulness of functional traits (i.e. estimated bite force) in reconstructing the
palaeo-ecology of shrews. We [1] investigate the relationships between mandible
morphology, bite force estimates, and the ecological context, and [2] assess the relevance of
estimated bite force as a palaeo-environmental indicator. To do so, we use geometric
morphometrics to quantify mandible shape diversity in shrews of the archaeological site El
Harhoura 2 (Rabat, Morocco), dated from the Late Pleistocene to the Holocene. To avoid
errors due to identification, morphological groups were used instead of taxa as units of
diversity. In order to explore how phenotypic traits are linked to their environment, mandible
shape and bite force estimates were compared with palaeo-environmental inferences for the
El Harhoura 2 site extracted from the literature. As hypothesized in previous studies,
morphological groups acted as phenotypic response units. Estimated bite force was related
to palaeo-environmental variation over the considered period, with a particular sensibility to
arid/humid environmental transitions. Functional characteristics of the morphological groups
pointed out a potential change in climate that had not been evidenced by previous palaeo-
environmental studies. The complementarity of morphological and functional indicators
allowed us to infer functional convergence and divergence among shrews. Our results
suggest that functional traits may be relevant indicators of changes in palaeo-environments.
However, they do not exhibit characteristic responses to particular climatic transitions. This
approach opens up new possibilities to explore the impact of environmental changes in
extinct organisms.
1. Introduction

Performance traits, provide a direct link between ecology, morphology, and fitness (Arnold,
1983; Wainwright, 1994). They reflect the ability of individuals to perform ecologically
relevant tasks (Irschick et al., 2008) and are subjected to selection (Irschick et al., 2008). In
animals, performance traits are dependent on skeletal structures, variation in muscular
anatomy, contractile physiology, and variation in biomechanical traits such as lever arms.
Consequently, the relationship between morphological and functional traits is complex
(Wainwright, 1994; Irschick et al., 2008). Different morphological traits can generate similar
functional outputs by redundancy (Alfaro, Bolnick, & Wainwright, 2005), which may lead to
functional convergence in organisms living in environments requiring similar performance
abilities (Wainwright, 2005; Young et al., 2007, 2010). Conversely, one morphological trait
can affect different performance traits through trade-offs and facilitation. Trade-offs occur
when there is a conflicting demand on a phenotypic trait through its differential implication in
several performance traits (Garland & Losos, 1994; Van Damme et al., 2003; Walker, 2007,
2010; Langerhans, 2009; Holzman et al., 2011; Vanhooydonck et al., 2011). Facilitation, on
the other hand, occurs when a similar demand is exerted on a phenotypic trait by several
different performance traits (Walker, 2007). Therefore, a complex "many-to-many"
morphology-performance relationship may exist (Bergmann & McElroy, 2014), implying that
morphological and functional traits may not respond similarly to changes in the environment.
Thus, functional traits may be relevant candidates to function as ecological indicators,
complementary to morphological traits.

Feeding is the function related to dietary ecology (Schwenk, 2000). Bite force is a
performance trait implicated in feeding. It is a performance trait directly linked to diet through
the mechanical demands imposed by variation in the mechanical properties of food items
(Anderson et al., 2008; Herrel et al., 2008; Dumont et al., 2009; Santana et al., 2010; Maestri
et al., 2016). Changes in diet associated with bite force are considered one of the main
drivers of diversification in mammals (Christiansen & Wroe, 2007; Monteiro & Nogueira,
2011). Thus, bite forces may be a powerful ecological proxy and have been widely used to
29 characterize dietary ecology in vertebrates (e.g. Freeman, 1979; Losos, 1992; Herrel et al., 2002; Huber, 2005; Kerr et al., 2017).

31 Studies on bite force in fossil taxa are common as they may provide insights into the
32 behavioural ecology of extinct taxa (e.g. Erickson et al., 1996; Therrien, 2005; Wroe et al.,
33 2005; Lappin et al., 2017; Rinderknecht et al., 2019). However, to our knowledge, estimates
34 of bite force have only rarely been applied to better understand palaeo-environmental
35 variation. Commonly, such inferences rely on community compositions and/or the
36 characterization of particular morphotypes (e.g. (Fernandez-Jalvo et al., 1998; Stoetzel et al.,
37 2011; Comay, Weissbrod & Dayan, 2019; Royer et al., 2020; López-García et al., 2021).
38 These approaches associate a particular biological feature to a specific environment.
39 However, an organism's dietary ecology can be more accurately inferred from bite force than
40 by categorical ecological classifications (Santana et al., 2010). Thus, bite force has a strong
41 potential to contribute meaningfully to understand paleo-environmental variation.

42 Functional aspects are also important determinants of morphological integration. This
43 concept defines the covariation patterns that exist between morphological traits (Badyaev &
44 Foresman, 2000, 2004; Young & Badyaev, 2006; Klingenberg, 2008; Klingenberg &
45 Marugán-Lobón, 2013). As a result of selection on performance, traits implicated in a same
46 function interact strongly, making them vary together (Klingenberg, 2008). This may be
47 particularly the case when they are linked through muscle insertions providing physical and
48 mechanical links between structures (Cheverud, Routman, & Irschick, 1997; Mezey,
49 Cheverud, & Wagner, 2000; Klingenberg, Mebus, & Auffray, 2003; Klingenberg, 2004). Thus,
50 changes in the strength of integration in the feeding apparatus could be related to changes in
51 bite force. Moreover, the less the modules in a structure like the mandible covary together
52 (i.e. the more modular they are), the more disparate are the organisms that can be produced
53 (Goswami et al., 2014). Modularity patterns may consequently also be related to bite force
54 disparity. Moreover, changes in modularity are known to respond to environmental stresses
55 (Badyaev & Foresman, 2000). In the light of these facts, modularity in the feeding apparatus
56 may also be a relevant indicator for understanding variation in palaeo-environments.
The Moroccan archaeological site El Harhoura 2 (EH2) is located in the Rabat-Témara region, on the North-Atlantic coast of Morocco. It is divided into 11 archeo-stratigraphical levels covering a time period from the Late Pleistocene to the Holocene. Phenotypic evolution related to environmental variation has been extensively studied at this site (Stoetzel, 2009; Stoetzel *et al.*, 2010, 2011, 2012b; Stoetzel, 2017; Stoetzel *et al.*, 2017). Small mammals such as rodents and shrews are a relevant model to consider palaeo-climatic variation (Valenzuela *et al.*, 2009; McGuire, 2010; Escudé *et al.*, 2013; Verde Arregoitia, Fisher & Schweizer, 2017). The ecological characteristics associated with species (as diet, habitat, soil) allowed palaeo-climatic inferences at EH2 (e.g. Stoetzel *et al.*, 2011, 2013, 2017). Among the remains recovered at EH2, shrew mandibles have been given peculiar attention (Cornette *et al.*, 2015a,c). Mandible shape is known for being highly plastic and often varies with environmental conditions in these animals (Badyaev & Foresman, 2000; Caumul & Polly, 2005; Cornette *et al.*, 2015a,c; Souquet *et al.*, 2019).

In the present study, we sought to understand the relationship between the morphology of the mandible in shrews, its modularity, and the estimated bite force over time. We further explored whether a functional trait (bite force) can provide relevant information for the understanding of palaeo-environments. Our approach is divided in three steps. 1) To evaluate global morphological diversity, shrews were divided into morphological groups (MGs). This characterization of biodiversity was chosen instead of species. MGs allow to consider variations at several levels of diversity (inter- and intra-specific), which makes it a relevant approach to characterize environment-related selection which acts at both levels (Erwin, 2000; Lande, 2009; Boutin & Lane, 2014; Hautmann, 2020). We expected identified MGs to be functional response units to changes in the environment, as has been suggested in previous studies (Read *et al.*, 2014; Khare *et al.*, 2017). 2) Bite force was estimated through the mechanical potential. Because performance and morphology are not strictly related, selection can apply differently at both levels of organization (Irschick *et al.*, 2008). Thus, quantifying variation in the mechanical potential may allow us to detect selection events undetectable through morphological proxies. 3) Finally, variation in modularity was
characterized. An increase in modularity in shrew mandibles may allow morphological diversification within a population utilizing the same resources irrespective of selection on function (Young et al., 2007, 2010) and as such be associated with variation in the environment.

89

## 2. Material and methods

### 2.1 *El Harhoura 2 cave and data collection*

El Harhoura 2 cave (33°57'08.9" N / 6°55'32.5" W) is an archaeological site located on the Moroccan Atlantic coast, a few km south of Rabat. The cave was occupied by Middle Stone Age, Later Stone Age, and Neolithic human populations and has yielded numerous archaeological material as well as the remains of large and small vertebrates (Nespoulet et al., 2008; El Hajraoui et al., 2012; Stoetzel et al., 2014). From top to bottom, the stratigraphy of EH2 cave is structured in 11 layers. Eight of these levels are well dated and were included in this study (L1, L2, L3, L4a, L5, L6, L7, L8) (Jacobs & Roberts, 2012; Ben Arous et al., 2019, 2020). Palaeo-environmental data were deduced from faunal communities (Stoetzel, 2009; Stoetzel et al., 2011, 2012b,a, 2014). An alternation of humid and arid periods has been documented and the landscape was mainly dominated by steppes with an increase in more wooded areas and water ponds during humid periods (Figure 1).

The material studied here is temporarily housed at the Muséum National d'Histoire Naturelle, Paris, France. It was sampled during the 2005-2009 excavation campaigns of the El Harhoura-Témara Archaeological Team (dir. R. Nespoulet and M.A. El Hajraoui). Four species of white-toothed shrews are represented: the material is largely dominated by *Crocidura russula* but *C. lusitania*, *C. tarfayensis* and *C. whitakeri* were also present (Cornette et al., 2015a,c). Among archaeological remains complete mandibles are rare and most of the material is fragmented. Nevertheless, mandible fragments also carry relevant morphological and palaeo-environmental information (Cornette et al., 2015a,c) and can be used to increase the sample. Here we used complete mandibles (Clpt) and three types of mandible fragments (A, B, C). Chosen fragments are those whose shape best enables
113

species to be distinguished (Cornette et al., 2015a). Fragmentation patterns are illustrated in **Figure 2**. Extant material (Act) of the same four species from the Rabat area was added to extend the time line to present day. The number of complete and fragmented mandibles studied is indicated in **Table 1**. Data acquisition is described in Cornette et al. (2015a).

114

115

116

117

118

## 2.2 Geometric morphometrics

119

Two dimensional geometric morphometric analyses were used. All mandibles could not be analyzed together because remains from various fragmentation patterns are not directly comparable. We performed separate shape analyses for each fragment type (Cplt, A, B and C) following the same protocol. Mandibles were analyzed through a landmark and sliding semi-landmark based approach allowing the description of the shape of biological relevant areas without anatomical landmarks (Bookstein, 1996; Zelditch, 2004; Gunz, Mitteroecker, & Bookstein, 2005; Cornette et al., 2013). Landmark locations for each fragment type are indicated in **Figure 2**. Semi-landmarks were slid to minimise the bending energy. Generalized Procrustes Analyses (GPA) were performed on each data set (Cplt, A, B and C) to remove effects of translation, rotation, and scale and make objects comparable (Rohlf & Slice, 1990). Resulting shape coordinates are the Procrustes residuals. These two last steps were performed using the "gpagen" function of the "Geomorph" package (Adams & Otárola-Castillo, 2013) in R version 3.6.3 (2020-02-29). To reduce data dimensionality, principal component analyses were performed on Procrustes residuals for each data set (Cplt, A, B and C) and we retained 95% of shape variation for the following analyses (Baylac & Frieß, 2005). All analyses were performed using R version 3.6.3 (2020-02-29).

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

## 2.3 Morphological groups

137

To assess shape diversity, mandibles were partitioned into morphological groups. First, complete mandibles (Cplt) were clustered based on their shape. To do so a morphological K-nearest neighbor method (KNN) was used. This method is adapted to small datasets, as is the case here, because it is a non-parametric classification. In this algorithm, each shape

138

139

140
object is assigned to its nearest neighbor cluster. The "clues" function of the "clues" package (Wang, Qiu & Zamar 2007) proposes an unsupervised KNN, meaning that the number of clusters is inferred from the data itself, favoring the most robust partitioning of the dataset. The robustness of the clusters is assessed by the Silhouette index (SI) (Kaufman & Rousseeuw, 1990) which measures the strength of the clusters. SI is comprised between -1 and 1. The more SI is close to 1 (i.e. SI>0), the more data points are correctly assigned to their clusters, and conversely the more SI is close to -1 (i.e. SI<0), the more data points are mis-assigned (Wang et al., 2007).

To visualize the morphological groups identified, the morphospace of complete mandibles (Cplt) was plotted using the three first axes of the principal component analysis computed on the Cplt data set. Deformations along axes were computed using the function "PlotRefToTarget" from the "geomorph" package (Adams & Otárola-Castillo, 2013). Thin-plate spline deformation grids representing differences between the extreme shapes of each axis and the global mean shape of Cplt mandibles were generated.

Next, we tested the robustness of the morphological groups identified on complete mandibles (Clpt) for each fragmentation pattern. Artificial A, B and C fragments were computed from complete mandibles (Clpt) and strengths of the clusters based on fragments were assessed using SI. This was achieved using the "get_Silhouette" function of the "clues" package (Wang et al., 2007).

Finally, the belonging of true A, B and C mandible fragments to the morphological groups was determined a posteriori using the KNN classification algorithm with the "knn" function of the "class" package (Venables & Ripley, 2002). Mean shapes of each morphological group were computed from complete mandibles (Clpt) using the functions "mshape" and "warpRefOutline" of the "Geomorph" package (Adams & Otárola-Castillo, 2013).

## 2.4 *Mechanical potential*

The mechanical potential as used in this study is the ratio of the muscle moment arm to the jaw outlever. As such it is dependent on the geometry of the skull and mandible and the
169 insertion of the masticatory muscles (Herrel et al., 2008; Chazeau et al., 2013; Manhães, 170 Nogueira, & Monteiro, 2017; Ginot et al., 2018, 2019). In particular, mandible shape is known 171 for being a good estimator of bite force (Brassard et al., 2020a,b). As a proxy for overall 172 mechanical potential, we choose the mechanical potential of the temporalis muscle. This is 173 one of the main muscles involved in bite force generation in shrews (Herrel et al., 2008; 174 Santana et al., 2010; Brassard et al., 2020a). It was estimated based on the moment arm of 175 the temporalis, computed from the mandible shape. The biomechanical model is presented 176 in Figure 3. It is defined as MP = A / B, where: A is the moment arm of the temporalis, B the 177 jaw out-lever, and MP the mechanical potential of the temporalis. MP was computed for 178 complete mandibles (Clpt) from Procrustes residuals.

179 However, MP does not reflect the phenotype in nature: because it is computed from 180 Procrustes residuals, it does not account for size variation in the data set, nor for 181 compensatory effects of musculature. Moreover, it only considers a single force - the 182 temporalis muscle - which is applied uniformly (Young et al., 2007). All these parameters are 183 important drivers of the mechanical potential, in particular size (Wroe et al., 2005; Freeman & 184 Lemen, 2008; Herrel et al., 2008; Chazeau et al., 2013; Manhães et al., 2017; Ginot et al., 185 2018, 2019; Brassard et al., 2020a). In order to increase the accuracy of our estimate of bite 186 force, we corrected MP to take size into account.

187 It is known that MP is linearly correlated to mandible/skull size in some mammals, 188 including shrews (Nogueira, Peracchi, & Monteiro, 2009; Cornette et al., 2015b; Manhães et 189 al., 2017; Ginot et al., 2018; Brassard et al., 2020a). As MP is computed from mandible 190 shape, this means that it has an allometric component, i.e. even if MP is size-free, a part of 191 the MP results from the influence of size. To confirm this assumption in our data set, we 192 tested it on complete mandibles (Cplt) by performing a regression of the log-transformed 193 centroid size on the log-transformed MP using the “lm” function of the “stats” package. 194 Centroid size (Csize) is a size estimator widely used in geometric morphometrics. It is 195 defined as the square root of squared distances of all landmarks of a mandible from its 196 centroid (Klingenberg, 2016).
197 Corrected MP (cMP) was subsequently expressed as: cMP = MP + f(size) (1), where MP
198 is the mechanical potential and f(size) the part of cMP due to size. It is also known that log₁₀
199 bite force is linearly correlated to size (Wroe et al., 2005; Chazeau et al., 2013; Manhães et
200 al., 2017; Ginot et al., 2019; Brassard et al., 2020a). Thus, for cMP to be a good estimator of
201 bite force, log-transformed cMP must also be linearly correlated to size. So, in (1) we have
202 f(size) = a(size + b (2)). The objective here is to determine a and b.

203 We assume that the allometric part of MP is proportional to the influence of size on MP.
204 Thus, a and b of (2) can be found by performing a linear regression of size on MP. However,
205 our data have a temporal component that we need to consider, otherwise, we may lose part
206 of this temporal information as MP and size are both related to time. We performed a
207 multivariate regression of size and time on MP. Log-transformed values of MP and Csize
208 were used. a and b of (2) were defined with the parameters of this regression: a as the
209 coefficient of log(Csize) and b as the intercept. We obtained: log(cMP) = log(MP) +
210 a*log(Csize) + b (3).
211

211 To evaluate whether cMP is a good estimator of bite force, we assessed its reliability on a
212 dataset of simulated data:

1. 1000 bite forces and associated sizes were simulated. Values were randomly generated according to a normal distribution using the "rnorm" function of the "stats" package, with for bite forces the constraint of a mean of 0.3 and a standard deviation of 0.15, and for size a mean of 3500 and a standard deviation of 500. Those constraints aimed to generate a dataset as close as possible to what is observed for the species present in our data set.

2. MP was computed using (3), where cMP was replaced by bite force values and with a
219 and b arbitrarily fixed (to represent the "real" relation between MP, bite force and size in the
220 simulated dataset). 30 couples of a and b were randomly generated using the "rnorm"
221 function of the "stats" package and tested. They were computed in order to be similar to the
222 a and b found previously. As the correlation between log₁₀ bite force and size is strictly
223 positive (Wroe et al., 2005; Chazeau et al., 2013; Manhães et al., 2017; Ginot et al., 2019;
224
Brassard et al., 2020a), a should be strictly positive and was generated with the constraint of a mean of 3 and a standard deviation of 3. b was generated with the constraint of a mean of 0 and a standard deviation of 5, to test positive and negative values. Tested couples of a and b are presented in Table S1. To simulate measurement error, Gaussian noise was added to computed MP (error tested at 0.01, 0.05 and 0.1).

3. cMP were calculated using (3), with a and b found previously.

4. For each pair of values in the simulated dataset, we tested whether the relation between the bite force of individuals was respected by cMPs (for example, if we have two individuals X and Y such that: bite force X > bite force Y, we must also have cMP X > cMP Y). Thus, variation in cMP reflects, in a relative way, variation in bite force. Mean, min and max scores of reliability among the 30 a and b couples were computed for the three measurement error estimates.

To obtain the cMP of fragmented mandibles (A, B and C) we predicted it from the cMP of complete ones (Clpt) according to the following protocol. First, to evaluate the reliability of the prediction, we performed covariation analyses between artificial fragments generated from the complete mandibles (Clpt) corresponding to fragmentation patterns (A, B and C) and cMP. We used two-block partial least squares analyses (PLS) which assesses the covariance between two sets of variables. This step was done using the "two.b.pls" function of the "geomorph" package (Adams & Otárola-Castillo, 2013), which performs a PLS adapted to shape data. Then, cMP of fragmented mandible was predicted using the "plsr" function of the "pls" package (Mevik & Wehrens, 2007). This function allows prediction based on the covariation between the two variables.

Shapes associated with the strongest and weakest cMP were computed. Those shapes were estimated based on the linear regression of cMP on shapes of complete mandibles (Cplt) using "procD.Im" of the "Geomorph" package (Adams & Otárola-Castillo, 2013). They were then computed with the "warpRefOutline" function of the same package.

Changes of cMP through time were investigated in two ways:
252

1) global cMP per layer. Differences between layers were tested through pairwise
comparisons testing using the "pairwise.t.test" function of the "stats" package which correct
for multiple testing.

253

254

255

2) cMP per morphological group per layer.

256

For both, standard deviation was computed on each layer as a measure of variance.

257

258

## 2.5 Modularity

259

As the abundance of complete mandibles (Clpt) per layer was not sufficient, the modularity analysis was performed based on the C fragmentation pattern on complete (Clpt) and C fragments. We divided the mandible into two modules *a priori* according to the hypothesized primary developmental modules defined in the literature (Cheverud *et al.*, 1997; Mezey *et al.*, 2000; Klingenberg *et al.*, 2003; Klingenberg, 2004) (**Figure 3**).

260

261

262

263

264

To quantify the modular structure of the mandible, we computed the Covariance Ratio (CR) per layer. CR compares the global covariation between hypothesized modules relatively to the covariation within those modules (Adams, 2016). The modularity hypothesis (independence of the hypothesized modules) is verified when CR<1. This measure is unaffected by sample size or the number of variables (Adams, 2016). It was performed using the "modularity.test" function of the "geomorph" package (Adams & Otárola-Castillo, 2013).

265

266

267

268

269

270

271

# 3. Results

272

## 3.1 Morphological groups

273

Three morphological groups were detected among complete mandibles (Clpt). We obtained an SI=0.15, meaning that the morphological partitions obtained with the clustering are correct. This partitioning was also valid for A, B, and C fragments (A fragments: SI=0.1; B fragments: SI=0.09; C fragments: SI=0.1).

274

275

276

277

A visualization of the morphological groups on the principal component analysis of complete mandibles (Cplt) and their deformations along major axes is presented in **Figure 4**.
Mean shapes of the morphological groups are presented in **Figure 5**. The first group

279
(illustrated in grey in **Figure 5**) has an elongated mandibular body, a more anteriorly inclined coronoid, and a more dorsally oriented condylar process compared to the average mandible shape. The second group (illustrated in red in **Figure 5**) displays a shape similar to the mean shape. The third group (illustrated in blue in **Figure 5**) has a short mandibular body, a coronoid process that is slightly more posteriorly inclined, and a condylar process that is more ventrally oriented.

286

## 3.2 Mechanical potential

Log-transformed MP revealed to be weakly and negatively correlated to log-transformed Csize (R²=0.014, P=0.0029**, 52 degrees of freedom). The multivariate regression (3) was also significant (R²=0.12, P=0.023*, 50 degrees of freedom). We obtained a=0.30 (coefficient of log-transformed Csize) and b=-3.47 (intercept). The reliability scores of cMP based on the simulated data are indicated in **Table 2**.

The PLS results showed that cMP is well related to A, B, and C fragments generated from complete mandibles (Clpt) (PLS: A fragments, r-pls=0.66, P<0.05; B fragments, r-pls=0.75, P<0.05; C fragments, r-pls=0.64, P<0.05), ensuring the reliability of predicted cMPs for fragments.

Mean shapes associated with the strongest and weakest cMP are presented in **Figure 6**. The shape associated with the weakest cMP (left in **Figure 6**) is characterized by an elongated mandibular body and a short mandibular branch compared to the global average shape of the mandible. On the contrary, the shape associated with the strongest cMP (right in **Figure 6**) displays a short mandibular body and a long mandibular branch.

Global changes in cMP through time are presented in **Figure 7A**. In L7, a period characterized by an open and arid environment, a great variability in cMPs co-exist compared to other layers. An important increase in average cMP is observed from L7 to L5 and then cMP decreases until present resulting in cMP values similar to those observed in L7. However, pairwise testing indicates that only the cMP of L5-L7 and Act-L5 are significantly different (P<0.05).

280

281

282

283

284

285

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307
Changes in cMP per morphological group through time are presented in **Figure 7B**. Not all morphological groups show the same variation in cMP over time, nor the same degree of variability. However, the three groups display a higher variability in L7. Overall, one of the morphological groups (indicated in blue in **Figure 7B**) has a higher cMP than the others. The two other groups (indicated in red and black in **Figure 7B**) display similar cMP values. In L5 the cMP of the three morphological groups converges toward high cMP values. Then, starting from L5, there is a drop of mandible shape diversity with the disappearance of one morphological group. In L4a and L3, the two remaining groups show a divergence in cMP towards respectively lower and higher values. In L2 a second diversity drop occurs with the disappearance of another group. From then onwards only one morphological group is consistently present, the two others show only few occurrences.

319

320 **3.3 Modularity**

Changes in CR over time are presented in **Figure 7C**. CR was not statistically significant in L2 and L4a. There are two main changes in CR values, first a slight increase in L7 (CR=0.86, P<0.05) followed by a strong increase from L5 to L3, with a peak in L3 (CR=0.95, P<0.05).

325

326 **4. Discussion**

In this study, we aimed [1] to investigate the relationship between shrew mandible morphology, modularity, and estimated biting forces (cMP) over time and [2] to evaluate whether a performance trait (bite force) could provide relevant information in the understanding of palaeo-environments. To do so, we explored mandible shape diversity in archaeological shrews, and computed the covariation ratio and estimated mechanical potential over time. First, we discuss the meaning of morphological groups and the benefits of this approach. Next, we focus on the relation between mandible shape and mechanical potential. Finally, we address the variation in these traits over time and compare them to

308

309

310

311

312

313

314

315

316

317

318

320

321

322

323

324

324

325

326

327

328

329

330

331

332

333

334
335

environmental inferences from the literature to assess their potential relevance for paleo-
environmental studies.

336

337

338

## 4.1 *Morphological groups*

339

Over the considered period, up to three morphological groups of mandibles were found among the four shrews species present at EH2 (Crocidura russula, C. whitakeri, C. tarfayensis and C. lusitania) (Cornette et al., 2015a,c). The assessed morphological diversity is consequently above the species level. Simplification of information is a risk of the morphological group approach (Read et al., 2014), but it provides access to a partition of variability that is potentially more informative ecologically than taxonomically when the studied biological object is sensitive to environmental change (Alperin, Cusminsky, & Bernasconi, 2011; Read et al., 2014; Khare et al., 2017), which is the case for shrew mandibles (Badyaev & Foresman, 2000; Young et al., 2010). The morphological groups reveal variations in the number of morphotypes independent of the number of species. For example, in the recent layers of EH2 (from L4a until today) only one or two morphological groups are present per layer whereas all four species are still present suggesting morphological convergence across these species (Cornette et al., 2015c). The morphological groups are characterized by different functional outputs and in their response to environmental changes. Except in L5, they display a different cMP and show different trends over time, especially from L4a and L3. To sum up, morphological groups are characterized by differences in morphology, functional output, and their response to environmental change. Thus, they represent morphological and functional response units to external variation.

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

Nevertheless, because we used only three of the five types of mandible fragments that were used in Cornette et al. (2015a, c), we reduced the initial sampling which can result in a potential loss of diversity. However, unused fragments were the least informative and reliable and could have introduced uncertainties into the results of this study which is why we chose to exclude them.

362
363 4.2 Mandible shape and mechanical potential

364 The relation between rostrum elongation and mechanical potential is intuitive as a longer
365 rostrum increases the jaw out-lever and consequently results in a lower mechanical potential.
366 This is illustrated by the negative correlation between MP and mandible size. This type of
367 relation between MP and mandible size has also been observed in other mammals (e.g.
368 Casanovas-Vilar & van Dam (2013) for squirrels, Nancy (1982) for felids). This implies that
369 the smallest specimens have a higher mechanical potential than the largest ones.
370 Mechanical potential and size are both important drivers of bite force. These results suggest
371 that the relative importance of these drivers varies between small and large specimens.
372 Conformations may be a more important driver of bite force than size in small specimens
373 compared to large ones. Nevertheless, the weakness of this correlation suggest that the
374 specimens used here are rather uniform in size.

375 Concerning the mandibular branch, our results are consistent with Young *et al.* (2007)
376 who found that a high mechanical potential was associated with a greater distance between
377 the condylar and the coronoid processes. This is related to muscles of the masticatory
378 apparatus: the coronoid process is the place of insertion of the temporalis muscle, the
379 condylar process of the external pterygoid and the angular process of the internal pterygoid
380 and the masseter. They all participate in generating bite force and impact bone shape in
381 shrews (Furió *et al.*, 2010; Cornette *et al.*, 2015c).

382 In the light of the previous, the mean shapes of morphological groups display
383 morphological features that can be related to their mechanical potential. The group
384 displaying the highest mechanical potential (in blue on **Figure 6** and **Figure 7**) is the one
385 with the shortest mandibular body and the closest condylar and coronoid processes, which
386 are morphological and functional particularities related to hard diet specialists in shrews
387 (Young *et al.*, 2007). Conversely the group with the weakest mechanical potential (in grey
388 **Figure 6** and in black on **Figure 7**) displays features characteristic of soft diet specialists in
389 shrews (Young *et al.*, 2007). The third group (in red on **Figure 6** and **Figure 7**) showing
390

average features likely regroups generalist shrews. Thus, morphological groups appear to
highlight ecological specialisations.

391

392

Modularity might be a key concept to understand the link between morphological and
functional variation. Three main decreases in modularity (i.e. increases in CR) are observed
in the sequence at EH2 in L7, L5, and L3. In L7, the three morphological groups display
unusual variability in cMP (it is important to note that this is the layer with the largest sample).
In L5, the three morphological groups show remarkably similar cMPs, meaning similar
functional outputs. Finally, in L3, the two morphological groups display highly divergent
functional outputs. Each decrease in modularity (i.e., increase in CR) in the mandible is
associated with an increase in the ability of a form to produce more diverse functional
outputs, allowing either divergence (as in L7 and L3) or convergence (as in L5) in the cMP.

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

These results are in contradiction with Young *et al.* (2010) and Young *et al.* (2007), who
found that extensive modularity allowed shrews with more diverse morphologies to produce a
similar functional output. However, Young *et al.* (2010) underlined that adaptive responses
are highly variable, even at a population level, which may explain the difference observed in
comparison to our results. Maybe the adaptive strategy here implies a different relation
between mechanical potential and the modularity of the mandible. A possibility might involve
variation in skull shape which was not studied here (Cornette *et al.*, 2015c).

408

409

4.3 *Mechanical potential as a palaeo-environmental indicator*

410

411

When considering both global and per morphological group changes in cMP, four
important functional variations were detected over time in L7, in L5, during L4a-L3 and in L2.
Three of these (in L7, L5 and L2) match transitions from humid to arid environments
(Stoetzel, 2009; El Hajraoui *et al.*, 2012).

412

413

414

415

416

417

In L7 we observed a high morphological diversity (with the presence of the three
morphological groups), the co-existence of highly diverse cMPs (Figure 7A and B) and a
decrease in modularity (Figure 7C). An increase in the diversity in cMP could be caused by a
release of selective pressures on this trait, allowing shrews with diverse abilities to survive.
This could be due to the availability of more diverse resources in the environment of L7 than in those of other layers. However, the environment in L7 is characterized by an open and dry environment of arid steppes, and seems to present less ecological diversity than some other layers (Stoetzel, 2009; El Hajraoui et al., 2012). Another hypothesis might be that a lack of resources drives character displacement in shrews with forms becoming highly specialized in the consumption of different resources. Moreover, the masseter and the medial pterygoid muscles also participate in generating bite force (Herrel et al., 2008; Santana et al., 2010; Brassard et al., 2020b,a). Yet, these muscles are implicated in the consumption of different types of resources, as in bats where the masseter allows the consumption of more soft resources (Santana et al., 2010). Moreover, functional variation is accompanied by a decrease in modularity. Variation in the degree of covariation between mandibular modules in shrews may be related to stressful environmental conditions (Badyaev & Foresman, 2000, p. 200). In the light of this fact, it is more likely that during the period covered by L7, shrews endured particularly stringent environmental conditions, with fewer and/or different available resources than before. As shrews are opportunistic, they might have switched their diet during this period inducing a release on the functioning of the temporalis muscle. However, as L7 is the layer with the largest sample this may bias our observations.

In L5, the cMP of two morphological groups greatly increases and the cMP of the three groups converges toward high values (Figure 7A and B). Pairwise testing indicates that global cMP is significantly different in L5 compared to L7 and present day (Figure 7A), and consequently supports the hypothesis of a cMP convergence towards high values in L5. This functional convergence is not associated with a convergence of the morphological groups, which is not surprising as those two types of convergence (functional and morphological) can be independent (Stayton, 2006). The group of soft diet shrews (in black on Figure 7B) is possibly subjected to a selection causing large sized specimens to produce the highest mechanical potential (as illustrated in Figure S1). Here we have a case of functional redundancy: three distinct morphologies producing similar functional outputs (Alfaro et al., 2005). This is observed in environmental conditions requiring similar performance abilities:
distinct morphologies are then able to adapt to similar functional demands (Wainwright, 2005; Young et al., 2007, 2010). This functional convergence could be explained by an increase in the selective pressures caused by fewer or different resources compared to previous layers, as the transition to L5 is towards a more arid environment (Stoetzel, 2009; El Hajraoui et al., 2012).

In L4a, the morphological group of soft diet shrews is not present anymore. It must have disappeared during the previous arid period during which soft resources may have been scarce. It is known that the insect cuticle becomes harder under dryer environments (Klocke & Schmitz, 2011). The cMP of the two remaining groups diverged distinctly. The supposed strong selection pressure(s) leading previously to the functional convergence in L5 must have eased. This is congruent with existing environmental inferences, as the environment in L4a is hypothesized to be very similar with conditions in L6 (Stoetzel, 2009; El Hajraoui et al., 2012). In L3, the divergence between cMP of the two groups increases. We can hypothesize that those two shrew morphotypes were subject to functional divergence. Each group of shrews may have specialized in the acquisition of a different food resource. The group displaying high cMP specializing in hard, large objects, requiring mechanical potential with an important contribution from the temporalis, and the other group specializing on softer, smaller items, which requires less contribution of the temporalis. This divergence might be caused by competition occurring among shrews. Such a competition has been suggested between *Crocidura russula* and other shrews species at EH2 at the exact same period (Cornette et al., 2015c), and shrew dietary specialization might be a response to competition (Smith & Remington, 1996). In L2, a second diversity drop is observed with the disappearance of the group of hard diet specialists. As in L5, this event occurred following an arid period during which less diverse resources may have been available. Only the group of generalist shrews is continuously present up to present.

Interestingly, extant shrews display a particularly weak cMP, without important environmental changes. Once again, this might be linked to a release of selective pressures on the temporalis resulting from a switch in diet (Santana et al., 2010). Another explanation
might be a selection towards weaker cMP. This is counter-intuitive, however, as selection for lower performance likely only occurs when the trait is energetically expensive to maintain or involved in trade-offs with other more relevant trait (Irschick et al., 2008). The energy previously allocated to the mechanical potential of the temporalis might have been reallocated to another performance trait under stronger selection in the novel environment.

## 5. Conclusion

To conclude, results of the present study illustrate the relevance of functional traits to detect palaeo-environmental transitions. Estimated bite forces showed variation in relation to palaeo-environmental changes over the considered period (Late Pleistocene - Holocene), with a particular sensibility to transition between arid and humid environments. The complementary nature of morphological and functional indicators allowed to infer and discuss the possible evolutionary and ecological processes involved. Functional traits have a great potential for refining paleo-environmental and palaeo-ecological inferences. Moreover, they appear to be relevant indicators palaeo-environmental transitions and offer a range of opportunities to explore the impact of environmental changes on extinct organisms.

## Acknowledgements

The authors would like to thank the Mission Archéologique El Harhoura-Témara, Ministère des Affaires Etrangères et Européennes, France, Ministère de la Culture, Morocco (dir. R. Nespoulet and M.A. El Hajraoui), the Institut National des Sciences de l'Archéologie et du Patrimoine (INSAP), Rabat, Morocco and the Institut Scientifique, Rabat, Morocco for helping to provide modern and fossil shrew specimens. We are very grateful to Peter Vogel for lending modern specimens of C. viaria and C. tarfayensis, and to Silvia Valenzuela-Lamas for lending specimens of C. russula from Spain. This research was funded by the Université Paris Descartes, the Ecole Doctorale FIRE – Programme Bettencourt. We also thank the two anonymous reviewers for constructive comments that helped to greatly improve the manuscript.
502

503 **Conflict of interest**

504 The authors declare that they have no conflict of interest.

505

506 **References**

507 **Adams DC. 2016.** Evaluating modularity in morphometric data: challenges with the RV coefficient and a new test measure (P Peres-Neto, Ed.). *Methods in Ecology and Evolution* **7**: 565–572.

510 **Adams DC & Otárola-Castillo E. 2013.** geomorph: an R package for the collection and analysis of geometric morphometric shape data (E Paradis, Ed.). *Methods in Ecology and Evolution* **4**: 393–399.

513 **Alfaro ME, Bolnick DI & Wainwright PC. 2005.** Evolutionary Consequences of Many-to-One Mapping of Jaw Morphology to Mechanics in Labrid Fishes. *The American Naturalist* **165**: E140–E154.

516 **Alperin MI, Cusminsky GC & Bernasconi E. 2011.** Benthic foraminiferal morphogroups on the argentine continental shelf. *The Journal of Foraminiferal Research* **41**: 155–166.

518 **Anderson RA, Mcbrayer LD & Herrel A. 2008.** Bite force in vertebrates: opportunities and caveats for use of a nonpareil whole-animal performance measure: bite force in vertebrates. *Biological Journal of the Linnean Society* **93**: 709–720.

521 **Arnold SJ. 1983.** Morphology, Performance and Fitness. *American Zoologist* **23**: 347–361.

522 **Badyaev AV & Foresman KR. 2000.** Extreme environmental change and evolution: stress-induced morphological variation is strongly concordant with patterns of evolutionary divergence in shrew mandibles. *Proceedings of the Royal Society of London. Series B: Biological Sciences* **267**: 371–377.

526 **Badyaev AV & Foresman KR. 2004.** Evolution of Morphological Integration. I. Functional Units Channel Stress-Induced Variation in Shrew Mandibles. *The American Naturalist* **163**: 868–879.

529 **Baylac M & Frieß M. 2005.** Fourier Descriptors, Procrustes Superimposition, and Data Dimensionality: An Example of Cranial Shape Analysis in Modern Human Populations. In: Slice DE, ed. *Modern Morphometrics in Physical Anthropology*. New York: Kluwer Academic Publishers-Plenum Publishers, 145–165.

533 **Ben Arous E, Falguères C, Tombret O, El Hajraoui MA & Nespoulet R. 2019.** Combined US-ESR dating of fossil teeth from El Harhoura 2 cave (Morocco): New data about the end of the MSA in Temara region. *Quaternary International*: S1040618218310668.

536 **Ben Arous E, Falguères C, Nespoulet R & El Hajraoui MA. 2020.** Review of chronological data from the Rabat-Temara caves (Morocco): Implications for understanding human occupation in north-west Africa during the Late Pleistocene. In: Leplongeon A, Goder-Goldberger M, Pleurdeau D, eds. *Not just a corridor. Human occupation of the Nile Valley and neighbouring regions between 75,000 and 15,000 years ago.*, 177–201.

542 **Bergmann PJ & McElroy EJ. 2014.** Many-to-Many Mapping of Phenotype to Performance:
543

An Extension of the F-Matrix for Studying Functional Complexity. *Evolutionary Biology*

544

41: 546–560.

545

Bookstein FL. 1996. Landmark Methods for Forms Without Landmarks: Localizing Group

546

Differences in Outline Shape. : 11.

547

Boutin S & Lane JE. 2014. Climate change and mammals: evolutionary versus plastic

548

responses. *Evolutionary Applications* 7: 29–41.

549

Brassard C, Merlin M, Guintard C, Monchâtre-Leroy E, Barrat J, Bausmayer N,

550

Bausmayer S, Bausmayer A, Beyer M, Varlet A, Houssin C, Callou C, Cornette R &

551

Herrel A. 2020a. Bite force and its relation to jaw shape in domestic dogs. *Journal of*

552

*Experimental Biology*: 47.

553

Brassard C, Merlin M, Monchâtre-Leroy E, Guintard C, Barrat J, Callou C, Cornette R

554

& Herrel A. 2020b. How Does Masticatory Muscle Architecture Covary with

555

Mandibular Shape in Domestic Dogs? *Evolutionary Biology* 47: 133–151.

556

Casanovas-Vilar I & van Dam J. 2013. Conservatism and Adaptability during Squirrel

557

Radiation: What Is Mandible Shape Telling Us? (A Goswami, Ed.). *PLoS ONE* 8:

558

e61298.

559

Caumul R & Polly D. 2005. Phylogenetic and environmental components of morphological

560

variation: skull, mandible and molar shape in marmots (Marmota, Rodentia). 59: 2460–

561

2472.

562

Chazeau C, Marchal J, Hackert R, Perret M & Herrel A. 2013. Proximate determinants of

563

bite force capacity in the mouse lemur: Bite force in mouse lemurs. *Journal of Zoology*

564

290: 42–48.

565

Cheverud JM, Routman EJ & Irschick DJ. 1997. Pleiotropic effects of individual gene loci

566

on mandibular morphology. *Evolution* 51: 2006–2016.

567

Christiansen P & Wroe S. 2007. Bite forces and evolutionary adaptations to feeding ecology

568

in Carnivores. *Ecology* 88: 347–358.

569

Comay O, Weissbrod L & Dayan T. 2019. Predictive modelling in paleoenvironmental

570

reconstruction: The micromammals of Manot Cave, Israel. *Journal of Human Evolution*:

571

102652.

572

Cornette R, Baylac M, Souter T & Herrel A. 2013. Does shape co-variation between the

573

skull and the mandible have functional consequences? A 3D approach for a 3D problem.

574

*Journal of Anatomy* 223: 329–336.

575

Cornette R, Herrel A, Stoetzel E, Moulin S, Hutterer R, Denys C & Baylac M. 2015a.

576

Specific information levels in relation to fragmentation patterns of shrew mandibles: do

577

fragments tell the same story? *Journal of Archaeological Science* 53: 323–330.

578

Cornette R, Tresset A, Houssin C, Pascal M & Herrel A. 2015b. Does bite force provide a

579

competitive advantage in shrews? The case of the greater white-toothed shrew: Bite

580

Force in Shrews. *Biological Journal of the Linnean Society* 114: 795–807.

581

Cornette R, Stoetzel E, Baylac M, Moulin S, Hutterer R, Nespoulet R, El Hajraoui MA,

582

Denys C & Herrel A. 2015c. Shrews of the genus Crocidura from El Harhoura 2

583

(Témara, Morocco): The contribution of broken specimens to the understanding of Late

584

Pleistocene–Holocene palaeoenvironments in North Africa. *Palaeogeography,*

585

*Palaeoclimatology, Palaeoecology* 436: 1–8.

586

Dumont ER, Herrel A, MedellÃ-n RA, Vargas-Contreras JA & Santana SE. 2009. Built
587 to bite: cranial design and function in the wrinkle-faced bat. *Journal of Zoology* **279**: 329–337.

589 **El Hajraoui MA**, **Nespoulet R**, **Debénath A** & **Dibble HL**. 2012. *Préhistoire de la région de Rabat-Témara*.

591 **Erickson GM**, **Kirk SDV**, **Su J**, **Levenston ME**, **Caler WE** & **Carter DR**. 1996. Bite-force estimation for Tyrannosaurus rex from tooth-marked bones. *Nature* **382**: 706–708.

593 **Erwin DH. 2000**. Macroevolution is more than repeated rounds of microevolution. *Evolution and Development* **2**: 78–84.

595 **Escudé É**, **Renvoisé É**, **Lhomme V** & **Montuire S. 2013**. Why all vole molars (Arvicolinae, Rodentia) are informative to be considered as proxy for Quaternary paleoenvironmental reconstructions. *Journal of Archaeological Science* **40**: 11–23.

598 **Fernandez-Jalvo Y**, **Denys C**, **Andrews P**, **Williams T**, **Dauphin Y** & **Humphrey L. 1998**. Taphonomy and palaeoecology of Olduvai Bed-I (Pleistocene, Tanzania). *Journal of Human Evolution* **34**: 137–172.

601 **Freeman PW. 1979**. Specialized Insectivory: Beetle-Eating and Moth-Eating Molossid Bats. *Journal of Mammalogy* **60**: 467–479.

603 **Freeman PW** & **Lemen CA. 2008**. A simple morphological predictor of bite force in rodents. *Journal of Zoology* **275**: 418–422.

605 **Furió M**, **Agustí J**, **Mouskhelishvili A**, **Sanisidro Ó** & **Santos-Cubedo A. 2010**. The paleobiology of the extinct venomous shrew *Beremendia* (Soricidae, Insectivora, Mammalia) in relation to the geology and paleoenvironment of Dmanisi (Early Pleistocene, Georgia). *Journal of Vertebrate Paleontology* **30**: 928–942.

609 **Garland T** & **Losos JB. 1994**. Ecological morphology of locomotor performance in squamate reptiles. In: Wainwright PC, Reilly SM, eds. *Ecological morphology: Integrative organismal biology.*, 240–302.

612 **Ginot S**, **Herrel A**, **Claude J** & **Hautier L. 2018**. Skull Size and Biomechanics are Good Estimators of *In Vivo* Bite Force in Murid Rodents: *in vivo* bite force estimation in murid rodents. *The Anatomical Record* **301**: 256–266.

615 **Ginot S**, **Herrel A**, **Claude J** & **Hautier L. 2019**. Morphometric models for estimating bite force in *Mus* and *Rattus*: mandible shape and size perform better than lever-arm ratios. *The Journal of Experimental Biology* **222**: jeb204867.

618 **Goswami A**, **Smaers JB**, **Soligo C** & **Polly PD. 2014**. The macroevolutionary consequences of phenotypic integration: from development to deep time. *Philosophical Transactions of the Royal Society B: Biological Sciences* **369**: 20130254–20130254.

621 **Gunz P**, **Mitteroecker P** & **Bookstein FL. 2005**. Semilandmarks in Three Dimensions. In: Slice DE, ed. *Modern Morphometrics in Physical Anthropology*. New York: Kluwer Academic Publishers-Plenum Publishers, 73–98.

624 **Hautmann M. 2020**. What is macroevolution? (A Smith, Ed.). *Palaeontology* **63**: 1–11.

625 **Herrel A**, **De Smet A**, **Aguirre LF** & **Aerts P. 2008**. Morphological and mechanical determinants of bite force in bats: do muscles matter? *Journal of Experimental Biology* **211**: 86–91.

628 **Herrel A**, **O'Reilly JC** & **Richmond AM. 2002**. Evolution of bite performance in turtles: Evolution of bite force in turtles. *Journal of Evolutionary Biology* **15**: 1083–1094.

630 **Holzman R**, **Collar DC**, **Mehta RS** & **Wainwright PC. 2011**. Functional Complexity Can
631 Mitigate Performance Trade-Offs. : 15.

632 **Huber DR. 2005.** Analysis of the bite force and mechanical design of the feeding mechanism
633 of the durophagous horn shark Heterodontus francisci. *Journal of Experimental Biology*
208: 3553–3571.

635 **Irschick DJ, Meyers JJ, Husak JF & Galliard JFL. 2008.** How does selection operate on
636 whole-organism functional performance capacities? A review and synthesis.
637 *Evolutionary Ecology Research* 10: 177–196.

638 **Jacobs Z, Roberts RG, Nespoulet R, El Hajraoui MA & Debénath A. 2012.** Single-grain
639 OSL chronologies for Middle Palaeolithic deposits at El Mnasra and El Harhoura 2,
640 Morocco: Implications for Late Pleistocene human-environment interactions along the
641 Atlantic coast of northwest Africa. *Journal of Human Evolution* 62: 377–394.

642 **Jacobs Z & Roberts RG. 2012.** Chapitre III. Datations par OSL avec la technique du grain
643 unique. In: *Préhistoire de la région de Rabat-Témara.*, 52–54.

644 **Janati-Idrissi N, Falgueres C, Nespoulet R, El Hajraoui MA, Debénath A, Bejjit L,  
645 Bahain JJ, Michel P, Garcia T, Boudad L, El Hammouti K & Oujaa A. 2012.**  
646 Datation par ESR-U/th combinées de dents fossiles des grottes d’El Mnasra et d’El  
647 Harhoura 2, région de Rabat-Temara. Implications chronologiques sur le peuplement du  
648 Maroc atlantique au Pléistocène supérieur et son. *Quaternaire* : 25–35.

649 **Kaufman L & Rousseeuw PJ. 1990.** *Finding groups in data. an introduction to cluster
650 analysis.*

651 **Kerr E, Cornette R, Rodrigues HG, Renaud S, Chevret P, Tresset A & Herrel A. 2017.**
652 Can functional traits help explain the coexistence of two species of Apodemus?
653 *Biological Journal of the Linnean Society:* 14.

654 **Khare N, Nigam R, Mayenkar DN & Saraswat R. 2017.** Cluster analysis of benthic
655 foraminiferal morpho-groups from the western margin of India reflects its depth
656 preference. *Continental Shelf Research* 151: 72–83.

657 **Klingenberg CP. 2004.** Integration and Modularity of Quantitative Trait Locus Effects on
658 Geometric Shape in the Mouse Mandible. *Genetics* 166: 1909–1921.

659 **Klingenberg CP. 2008.** Morphological Integration and Developmental Modularity. *Annual
660 Review of Ecology, Evolution, and Systematics* 39: 115–132.

661 **Klingenberg CP. 2016.** Size, shape, and form: concepts of allometry in geometric
662 morphometrics. *Development Genes and Evolution* 226: 113–137.

663 **Klingenberg CP & Marugán-Lobón J. 2013.** Evolutionary Covariation in Geometric
664 Morphometric Data: Analyzing Integration, Modularity, and Allometry in a Phylogenetic
665 Context. *Systematic Biology* 62: 591–610.

666 **Klingenberg CP, Mebus K & Auffray JC. 2003.** Developmental integration in a complex
667 morphological structure: how distinct are the modules in the mouse mandible? *Evolution
668 and Development* 5: 522–531.

669 **Klocke D & Schmitz H. 2011.** Water as a major modulator of the mechanical properties of
670 insect cuticle. *Acta Biomaterialia* 7: 2935–2942.

671 **Lande R. 2009.** Adaptation to an extraordinary environment by evolution of phenotypic
672 plasticity and genetic assimilation. *Journal of Evolutionary Biology* 22: 1435–1446.

673 **Langerhans RB. 2009.** Trade-off between steady and unsteady swimming underlies predator-
674 driven divergence in *Gambusia affinis*. *Journal of Evolutionary Biology* 22: 1057–1075.
675 Lappin AK, Wilcox SC, Moriarty DJ, Stoeppler SAR, Evans SE & Jones MEH. 2017.

676 Bite force in the horned frog (Ceratophrys cranwelli) with implications for extinct giant

677 frogs. *Scientific Reports* 7: 11963.

678 López-García JM, Cuenca-Bescós G, Galindo-Pellicena MÁ, Luzi E, Berto C, Lebreton

679 L & Desclaux E. 2021. Rodents as indicators of the climatic conditions during the

680 Middle Pleistocene in the southwestern Mediterranean region: implications for the

681 environment in which hominins lived. *Journal of Human Evolution* 150: 102911.

682 Losos JB. 1992. The evolution of convergent structure in caribbean anolis communities.

683 *Systematic Biology* 41: 18.

684 Maestri R, Patterson BD, Fornel R, Monteiro LR & de Freitas TRO. 2016. Diet, bite

685 force and skull morphology in the generalist rodent morphotype. *Journal of Evolutionary

686 *Biology* 29: 2191-2204.

687 Manhães IA, Nogueira MR & Monteiro LR. 2017. Bite force and evolutionary studies in

688 phyllostomid bats: a meta-analysis and validation. *Journal of Zoology* **302**: 288-297.

689 McGuire JL. 2010. Geometric morphometrics of vole (Microtus californicus) dentition as a

690 new paleoclimate proxy: Shape change along geographic and climatic clines. *Quaternary*

691 *International* **212**: 198-205.

692 Mevik BH & Wehrens R. 2007. The *pls* Package: Principal Component and Partial Least

693 Squares Regression in *R. Journal of Statistical Software* **18**.

694 Mezey JG, Cheverud JM & Wagner GP. 2000. Is the Genotype-Phenotype Map Modular?:

695 A Statistical Approach Using Mouse Quantitative Trait Loci Data. : 7.

696 Monteiro LR & Nogueira MR. 2011. Evolutionary patterns and processes in the radiation of

697 phyllostomid bats. *BMC Evolutionary Biology* **11**: 137.

698 Nancy A. N. 1982. *The Big Cats: The Paintings of Guy Coheleach*.

699 Nespoulet R, El Hajraoui MA, Amani F, Ben Ncer A, Debénath A, El Idrissi A,

700 Lacombe JP, Michel P, Oujaa A & Stoetzel E. 2008. Palaeolithic and Neolithic

701 Occupations in the Témara Region (Rabat, Morocco): Recent Data on Hominin Contexts

702 and Behavior. *African Archaeological Review* **25**: 21-39.

703 Nespoulet R & El Hajraoui MA. 2012. Excavation report.

704 Nogueira MR, Peracchi AL & Monteiro LR. 2009. Morphological correlates of bite force

705 and diet in the skull and mandible of phyllostomid bats. *Functional Ecology* **23**: 715-723.

706 Read CF, Duncan DH, Vesk PA & Elith J. 2014. Biocrust morphogroups provide an

707 effective and rapid assessment tool for drylands (S Wan, Ed.). *Journal of Applied*

708 *Ecology* **51**: 1740-1749.

709 Rinderknecht A, Jones WW, Araújo N, Grinspan G & Blanco RE. 2019. Bite force and

710 body mass of the fossil rodent *Telecomys giganteus* (Caviomorpha, Dinomyidae).

711 *Historical Biology* **31**: 644-652.

712 Rohlf FJ & Slice D. 1990. Extensions of the Procrustes Method for the Optimal

713 Superimposition of Landmarks. *Systematic Zoology* **39**: 40.

714 Royer A, Garcia Yelo BA, Laffont R & Hernandez Fernandez M. 2020. New bioclimatic

715 models for the quaternary palaearctic based on insectivore and rodent communities.

716 *Palaeogeography, Palaeoclimatology, Palaeoecology* **560**: 18.

717 Santana SE, Dumont ER & Davis JL. 2010. Mechanics of bite force production and its

718 relationship to diet in bats: Bite force and diet in bats. *Functional Ecology* **24**: 776-784.
719 **Schwenk K. 2000.** *Form, Function, and Evolution in Tetrapod Vertebrates.*

720 **Smith HR & Remington CL. 1996.** Food Specificity in Interspecies Competition.

721 *BioScience* **46**: 436–447.

722 **Souquet L, Chevret P, Ganem G, Auffray JC, Ledevin R, Agret S, Hautier L & Renaud**

723 S. **2019.** Back to the wild: does feralization affect the mandible of non-commensal house mice (*Mus musculus domesticus*)? *Biological Journal of the Linnean Society* **126**: 471–486.

726 **Stayton CT. 2006.** Testing hypotheses of convergence with multivariate data: morphological and functional convergence among herbivorous lizards. *Evolution* **60**: 824–841.

728 **Stoetzel E. 2009.** Les microvertébrés du site d'occupation humaine d'El Harhoura 2 (Pleistocene supérieur - Holocène, Maroc) : systématique, évolution, taphonomie et paléoécologie. : 437.

731 **Stoetzel E, Bailon S, Nespoulet R, El Hajraoui MA & Denys C. 2010.** Pleistocene and holocene small vertebrates of El Harhoura 2 cave (Rabat-Témara, Morocco): an annotated preliminary taxonomic list. *Historical Biology* **22**: 303–319.

734 **Stoetzel E, Marion L, Nespoulet R, El Hajraoui MA & Denys C. 2011.** Taphonomy and palaeoecology of the late Pleistocene to middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Témara, Morocco). *Journal of Human Evolution* **60**: 1–33.

737 **Stoetzel E, Bougariane B, Campmas E, Ouchaou B & Michel P. 2012a.** Chapitre V. Faunes et paléoenvironnements. In: *Préhistoire de la région de Rabat-Témara*, 35–51.

739 **Stoetzel E, Denys C, Bailon S, El Hajraoui MA & Nespoulet R. 2012b.** Taphonomic Analysis of Amphibian and Squamate Remains from El Harhoura 2 (Rabat-Témara, Morocco): Contributions to Palaeoecological and Archaeological Interpretations: Taphonomic Study of Amphibian and Squamate Fossil Remains. *International Journal of Osteoarchaeology* **22**: 616–635.

744 **Stoetzel E, Denys C, Michaux J & Renaud S. 2013.** *Mus* in Morocco: a Quaternary sequence of intraspecific evolution: Quaternary Evolution of *Mus* in Morocco. *Biological Journal of the Linnean Society* **109**: 599–621.

747 **Stoetzel E, Campmas E, Michel P, Bougariane B, Ouchaou B, Amani F, El Hajraoui MA & Nespoulet R. 2014.** Context of modern human occupations in North Africa: Contribution of the Témara caves data. *Quaternary International* **320**: 143–161.

750 **Stoetzel E, Cornette R, Lalis A, Nicolas V, Cucchi T & Denys C. 2017.** 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. *Quaternary Science Reviews* **164**: 199–216.

754 **Stoetzel E. 2017.** Adaptations and Dispersals of Anatomically Modern Humans in the Changing Environments of North Africa: the Contribution of Microvertebrates. *African Archaeological Review* **34**: 453–468.

757 **Therrien F. 2005.** Feeding behaviour and bite force of sabretoothed predators. *Zoological Journal of the Linnean Society* **145**: 393–426.

759 **Valenzuela S, Poitevin F, Cornette R, Bournery A, Nadal J & Vigne JD. 2009.** Evolving ecosystems: ecological data from an Iron Age small mammal accumulation at Alorda Park (Catalonia, Spain). *Journal of Archaeological Science* **36**: 1248–1255.

762 **Van Damme R, Vanhooydonck B, Aerts P & De Vree F. 2003.** Evolution of lizard
763 locomotion: Context and constraint. In: *Vertebrate Biomechanics and Evolution*, 267–

764 282.

765 **Vanhooydonck B, Boistel R, Fernandez V & Herrel A. 2011.** Push and bite: trade-offs

766 between burrowing and biting in a burrowing skink (Acontias percivali): performance

767 trade-offs in a burrowing skink. *Biological Journal of the Linnean Society* **102**: 91–99.

768 **Venables WN & Ripley BD. 2002.** Modern Applied Statistics with S. : 504.

769 **Verde Arregoitia LD, Fisher DO & Schweizer M. 2017.** Morphology captures diet and

770 locomotor types in rodents. *Royal Society Open Science* **4**: 160957.

771 **Wainwright PC. 1994.** Functional morphology as a tool in ecological research. In:

772 *Ecological Morphology*. (eds P.C. Wainwright & S.M. Reilly)., 42–59.

773 **Wainwright PC. 2005.** Many-to-One Mapping of Form to Function: A General Principle in

774 Organismal Design? *Integrative and Comparative Biology* **45**: 256–262.

775 **Walker JA. 2007.** A General Model of Functional Constraints on Phenotypic Evolution. *The*

776 *American Naturalist* **170**: 681–689.

777 **Walker JA. 2010.** An Integrative Model of Evolutionary Covariance: A Symposium on Body

778 Shape in Fishes. *Integrative and Comparative Biology* **50**: 1051–1056.

779 **Wang X, Qiu W & Zamar RH. 2007.** CLUES: A non-parametric clustering method based

780 on local shrinking. *Computational Statistics & Data Analysis* **52**: 286–298.

781 **Wroe S, McHenry C & Thomason J. 2005.** Bite club: comparative bite force in big biting

782 mammals and the prediction of predatory behaviour in fossil taxa. *Proceedings of the*

783 *Royal Society B: Biological Sciences* **272**: 619–625.

784 **Young RL & Badyaev AV. 2006.** Evolutionary persistence of phenotypic integration:

785 influence of developmental and functional relationships on complex trait evolution.

786 *Evolution* **60**: 1291–1299.

787 **Young RL, Haselkorn TS & Badyaev AV. 2007.** Functional equivalence of morphologies

788 enables morphological and ecological diversity. *Evolution* **61**: 2480–2492.

789 **Young RL, Sweeney MJ & Badyaev AV. 2010.** Morphological diversity and ecological

790 similarity: versatility of muscular and skeletal morphologies enables ecological

791 convergence in shrews: *Morphological versatility and evolution. Functional Ecology* **24**:

792 556–565.

793 **Zelditch M (Ed.). 2004.** *Geometric morphometrics for biologists: a primer*. Amsterdam ;

794 Boston: Elsevier Academic Press.

795

796
797 **Figure captions**

798 **Figure 1:** Stratigraphy, datation and characteristic environments of the eight layers of EH2
799 used in this study (Stoetzel, 2009; Jacobs et al., 2012; Janati-Idrissi et al., 2012; Nespoulet &
800 El Hajraoui, 2012; Stoetzel et al., 2012a; Ben Arous et al., 2019, 2020). Three datations
801 methods were used: AMS-14C (based on organic remains), US-ESR (Combined uranium
802 series and Electron Spin Resonance) and OSL (Optical Stimulated Luminescence). Light
803 blue layers: humid and open environment; dark blue layers: humid and close environment;
804 light yellow layers: arid and open environment; dark yellow layers: arid and semi-open
805 environment. Unused layers in this study are in grey.

806 **Figure 2:** Types of fragments and their landmark (red points) and semi-landmark (blue
807 points) locations. Clpt: complete mandibles; A: fragmented mandibles of type A; B:
808 fragmented mandibles of type B; C: fragmented mandibles of type C.

809 **Figure 3:** Biomechanical model used to estimate the mechanical potential of the temporalis
810 muscle (MP) from mandible shape. V is the vector that starts at mid distance between
811 landmarks 3 and 5 and has for direction the intersection between the line passing through
812 landmarks 3 and 4 and the parallel of the line passing through landmarks 2 and 3 passing
813 through 5. A is the moment arm if the temporalis (the distance between the landmark 2 and
814 the vector). B is the jaw outlever (the distance between landmarks 1 and 2). M1 and M2 are
815 the definitions of the two hypothesized modules of the mandible from the literature (Cheverud
816 et al., 1997; Mezey et al., 2000; Klingenberg et al., 2003; Klingenberg, 2004) and landmarks
817 implicated (C fragmentation pattern). M1: ascending ramus module; M2: alveolar region
818 module.

819 **Figure 4:** Principal component analysis on complete mandibles (Cplt).

820 **Figure 5:** Mean shapes of the three morphological groups (in grey, red and blue) compared
821 to the global mean shape (in black) of complete mandibles (Clpt).

822 **Figure 6:** Mean shapes associated with the strongest (right) and weakest (left) cMP (in
823 yellow) compared to the global mean shape (in black) of complete mandibles (Clpt). Strength
824 of the cMP is symbolized by the size of the shrew drawing under the models.
**Figure 7:** Global changes in cMP (A), cMP per morphological groups (B) and CR (C) of the shrews of EH2, from L8 until current day. Environmental conditions (Stoetzel, 2009; El Hajraoui *et al.*, 2012) are indicated by background colors. For cMP, standard deviation is indicated for each point. A: layers displaying significant different cMP are indicated in red. B: Black is the first morphological group; red is the second morphological group; blue is the third morphological group. C: Dotted lines indicate missing values.
831 **Table 1:** Abundance of remains studied (complete and fragmented mandibles).

<table><thead><tr><th>Layer</th><th>Clpt</th><th>A</th><th>B</th><th>C</th><th>TOTAL</th></tr></thead><tbody><tr><td>Act</td><td>16</td><td>-</td><td>-</td><td>-</td><td>16</td></tr><tr><td>L1</td><td>2</td><td>1</td><td>-</td><td>2</td><td>5</td></tr><tr><td>L2</td><td>2</td><td>4</td><td>1</td><td>2</td><td>9</td></tr><tr><td>L3</td><td>3</td><td>1</td><td>1</td><td>1</td><td>6</td></tr><tr><td>L4a</td><td>2</td><td>1</td><td>-</td><td>1</td><td>4</td></tr><tr><td>L5</td><td>10</td><td>-</td><td>2</td><td>3</td><td>15</td></tr><tr><td>L6</td><td>10</td><td>2</td><td>4</td><td>9</td><td>25</td></tr><tr><td>L7</td><td>9</td><td>7</td><td>13</td><td>10</td><td>39</td></tr><tr><td>L8</td><td>-</td><td>6</td><td>10</td><td>11</td><td>27</td></tr><tr><td>TOTAL</td><td></td><td></td><td></td><td></td><td>146</td></tr></tbody></table>

832

833
834

<table><caption>**Table 2:** Reliability scores of cMP for the three tested measure errors.</caption><thead><tr><th>Measure error</th><th>0.01</th><th>0.05</th><th>0.1</th></tr></thead><tbody><tr><th>Mean</th><td>82%</td><td>82%</td><td>81%</td></tr><tr><th>Max</th><td>98%</td><td>96%</td><td>93%</td></tr><tr><th>Min</th><td>66%</td><td>66%</td><td>66%</td></tr></tbody></table>

835

836
838 **Table S1:** Couples of *a* and *b* tested on the simulated data set of bite forces and cMPs to
839 assess the reliability of cMP as an estimator of bite force.

<table><thead><tr><th>a</th><th>b</th></tr></thead><tbody><tr><td>3.64</td><td>0.87</td></tr><tr><td>6.97</td><td>9.11</td></tr><tr><td>2.30</td><td>4.38</td></tr><tr><td>1.60</td><td>2.63</td></tr><tr><td>5.74</td><td>-2.40</td></tr><tr><td>8.56</td><td>-7.38</td></tr><tr><td>0.89</td><td>-1.81</td></tr><tr><td>7.09</td><td>5.49</td></tr><tr><td>1.89</td><td>-1.52</td></tr><tr><td>1.35</td><td>5.30</td></tr><tr><td>0.42</td><td>2.27</td></tr><tr><td>4.44</td><td>1.08</td></tr><tr><td>2.21</td><td>4.32</td></tr><tr><td>2.01</td><td>-6.27</td></tr><tr><td>5.74</td><td>-1.78</td></tr><tr><td>0.92</td><td>3.07</td></tr><tr><td>4.55</td><td>-5.31</td></tr><tr><td>3.34</td><td>0.82</td></tr><tr><td>5.90</td><td>0.86</td></tr><tr><td>7.90</td><td>4.29</td></tr><tr><td>0.40</td><td>-3.71</td></tr><tr><td>1.64</td><td>-7.20</td></tr><tr><td>1.03</td><td>9.24</td></tr></tbody></table>
<table><tr><td>0.38</td><td>-2.87</td></tr><tr><td>2.82</td><td>-4.73</td></tr><tr><td>4.60</td><td>-3.80</td></tr><tr><td>1.64</td><td>1.13</td></tr><tr><td>4.73</td><td>2.03</td></tr><tr><td>5.78</td><td>1.55</td></tr><tr><td>5.58</td><td>-0.37</td></tr></table>
**Figure S1:** Global changes in Csize (A) and Csize per morphological groups (B) of the shrews of EH2, from L8 until current day. Environmental conditions (Stoetzel, 2009; El Hajraoui *et al.*, 2012) are indicated by background colors. Standard deviation is indicated for each point. B: Black is the first morphological group; red is the second morphological group; blue is the third morphological group.

**Figure S2:** Global changes in the moment arm of the temporalis (A on Figure 3) (A) and lever arm of the temporalis per morphological groups (B) of the shrews of EH2, from L8 until current day. Environmental conditions (Stoetzel, 2009; El Hajraoui *et al.*, 2012) are indicated by background colors. Standard deviation is indicated for each point. B: Black is the first morphological group; red is the second morphological group; blue is the third morphological group.