Trends in functional composition of small mammal communities across millennial time scales

Rich fossil deposits of the late Quaternary help us understand responses of biodiversity to global change and thus predict the future of ecosystems. Studies from the late Quaternary, however, are often limited taxonomically, geographically (often one site), and by their use of largely taxon-based metrics that do not inform about ecosystem-level consequences of biodiversity change. Here, we compare change in functional composition of small mammal communities at El Mirón Cave (Spain) and Samwell Cave (California, USA) across the last 22 000 years, and examine their relationships with climate and vegetation. We find opposing temporal trends between the two locations. European small mammal communities occupied increasingly greater trait space, driven by increases in arboreal granivory and frugivory as ground-dwelling herbivory declined toward the late Holocene. North American communities occupied smaller trait space as ground-level foraging and insectivory increased and arboreal herbivory and mean body mass declined. Our results point to the importance of the interaction between climate change and vegetation shifts for explaining changes in small mammal functional diversity, through their synergistic impacts on individual traits. Specifically, increasing temperature across both continents likely led to increases in nocturnal activity and declines in assemblage mean body mass, while transition to mixed forest (El Mirón) or open woodland (Samwell) resulted in increasing structural complexity of vegetation that potentially supported more diverse community-level dietary characteristics. These results demonstrate the ability of a trait-based approach to identify how environmental variables correlate with changes in community functional composition through time and gain insight into the potential consequences of environmental change for ecosystem functioning.


Introduction
The Earth is undergoing a sixth mass extinction event (Barnosky et al. 2011b, Ceballos et al. 2020) and redistribution of species (Pecl et al. 2017) as a result of anthropogenic change (Díaz et al. 2019). This predicted loss and ongoing redistribution of life on Earth will have staggering consequences for communities, ecosystems and humans (Pecl et al. 2017, IPBES 2019. Responses to global change at the community and ecosystem levels, however, remain poorly understood, hindering sound forecasting of future states of biodiversity and selecting appropriate conservation actions. Our understanding is limited in part because typical scales of community ecology research that these predictions rely on are often insufficient to understand the long-term consequences of species redistribution and extinction (Estes et al. 2018). To elucidate the effects of global change on ecological systems, investigations must be carried out at scales consistent with those of climate change, including time scales much longer than human lifetime (Barnosky et al. 2011b).
The late Quaternary (i.e. the late Pleistocene and Holocene, from the Last Glacial Maximum [LGM] ca 22 000 years before present [22 kybp] to today) provides an excellent opportunity to explore responses of biota to global change and help predict the future of ecosystems in the face of increasing human pressure (Gill et al. 2009, Blois et al. 2010, Jackson and Blois 2015, Lyons et al. 2016, Davis 2017, Smith et al. 2018, Tóth et al. 2019. The Pleistocene-Holocene transition includes a time of rapid climate change Overpeck 2000, Broecker et al. 2010) with comparable magnitude to that of predicted anthropogenic climate change over the next few centuries (Jackson and Overpeck 2000). Furthermore, the late Quaternary fossil record is exemplary in its completeness, quality and quantity and is dated with relatively high accuracy and precision (Jackson and Blois 2015), allowing for fine resolution assessment of ecological dynamics.
Fossil pollen and large-bodied mammals are predominantly represented in studies from the late Quaternary, but they may not by themselves provide a complete understanding of the impacts of global change on ecological dynamics and their consequences for ecosystem functioning. Plants are mostly restricted to the foundation of food webs, which limits the dimensionality of potential ecological functions to those provided by autotrophs (Pigot et al. 2020) and might underestimate the full range and significance of changes in communities to ecosystems. Likewise, fossil deposits of large-bodied mammals may not always reflect the composition of entire mammalian assemblages (Behrensmeyer et al. 1979) and might thus miss crucial functional characteristics of assemblages. Small mammals fossil deposits, on the other hand, preserve animals that constitute an interacting community (McGill et al. 2005) that is spatially and temporally concordant (Porder et al. 2003, Pardi andGraham 2018). Small mammals additionally display high climate and vegetation fidelity (Hadly 1996, Moritz et al. 2008, Rowe et al. 2011, Rowe and Terry 2014, respond rapidly to environmental and biotic change (Rowe and Terry 2014), and provide both direct (Brown and Heske 1990, Keesing 2000, Bricker et al. 2010) and indirect (Root-Bernstein et al. 2013) ecological functions, thus having the potential to provide a clear link between environmental changes and changes in ecosystem functioning.
Previous work on small mammals has identified correlations between their diversity dynamics and the late Quaternary environmental changes (Graham 1976, Grayson 2000, Carrasco et al. 2009, Blois et al. 2010, Terry et al. 2011), but few have made explicit links to ecosystem functions. By being descriptors of species' niches and their functional roles, functional and life history traits are emerging as a primary way to gain insight into ecosystem-level consequences of biodiversity change. Furthermore, by being directly linked to environmental gradients (McGill et al. 2006), traits also provide a better understanding of ecosystem functionality in the face of environmental change (Pigot et al. 2020). Increasing availability of ecological trait data now allows a more mechanistic, taxon-free evaluation of biodiversity dynamics (McGill et al. 2006, Lawing et al. 2012, Polly and Head 2015, Jarzyna and Jetz 2016. To date, only a handful of studies have applied a trait-based approach to the late Quaternary and, as with most research on biotic change during this time period, they are limited to North America (Terry et al. 2011, Davis 2017. Given that changes in climate and vegetation were not temporally synchronous or of equal magnitude across continents (Barnosky et al. 2004, Smith et al. 2018), a geographic bias toward North America hinders attaining a more complete understanding of the ecosystem-level consequences of biodiversity change by limiting our knowledge of biotic response to just one of myriad contexts.
Here, we investigate changes in functional composition of small mammal communities at two sites: El Mirón Cave in northern Spain (Cuenca-Bescós et al. 2009) and Samwell Cave in northern California (Blois et al. 2010). El Mirón and Samwell Caves underwent significant climatic and vegetation changes over the past 20 000 years and are representative of the larger regions in which they are located. Both sites capture at least 15 000 years of biotic change in stratigraphically distinguishable and radiocarbon dated small mammal communities, and as such provide some of the most comprehensive and highest resolution record of small mammal communities across the late Quaternary. To gain a comprehensive and multi-dimensional understanding of changes in mammalian trait space (i.e. multi-dimensional space wherein species are placed according to their combinations of traits; Mouillot et al. 2021), we use a suite of multivariate and univariate functional and trait diversity metrics. We further link changes in each of the functional and trait diversity metrics to region-specific environmental (i.e. climate and vegetation) conditions and draw conclusions regarding the ecosystem-level consequences of observed changes in small mammal diversity. We predict that widespread environmental change led to significant changes in functional composition of small mammal assemblages by promoting select trait characteristics.

Chronology and abundance data
To compare change in the functional composition of small mammal communities in Europe and North America, we used abundance and chronology data for two sites: El Mirón Cave (EMC) in Spain (Cuenca-Bescós et al. 2009) and Samwell Cave (SC) in California (Blois et al. 2010). Fossil data collection at both sites was done using screenwashing. SC data were downloaded from the Neotoma Paleoecology Database (Williams et al. 2018) using the R package 'neotoma' ver. 1.7.4 (Goring et al. 2015) and contain 14 layers with associated calibrated radiocarbon dates, ranging from 0.75 to 17.7 kybp. EMC data contains 32 layers with associated radiocarbon dates from Cuenca-Bescós et al. (2009) and additional dates from Straus and González Morales (2003 and Straus et al. (2015). We removed the two oldest dated layers due to large time gaps between them and the next youngest dates, potentially obscuring our ability to detect trends, and the most recent layer due to an insufficient number of specimens (two Talpa europea and two Lagomorpha indet). We further excluded layer 115 from the analysis due to its imprecise age estimate (Straus and González Morales 2018). Ultimately, 28 dated layers at EMC remained available for analysis. Each radiocarbon date at EMC was calibrated using the 'BchronCalibrate' function in the package 'Bchron' ver. 4.3.0 using the 'intcal13' calibration curve (Haslett and Parnell 2008). A distribution of 20 000 ages for each layer at EMC was produced using the function 'sampleAges', and the median value was extracted from this distribution. For layers with multiple ages associated with them, we used the mean of the median ages.
As is typical for paleontological studies of small mammals, we excluded bats from the analysis due to the taphonomic differences in preservation between bats and nonvolant terrestrial mammals (Andrews 1990), and removed all occurrences not resolved to at least the genus level. Most occurrences at EMC were resolved to the species level. The taxonomic list includes occurrences for Apodemus gr. sylvaticus-flavicollis and Sorex gr. coronatus-araneus, and for these occurrences, we used mean trait values of the two species. Nine individuals were resolved only to the order Lagomorpha and were removed from the species list. Pliomys lenki, which accounted for 33 occurrences (22 of which were found in the earliest, discarded, layer), was removed from the dataset because it is extinct and thus has no available trait data. At SC, most occurrences (3228 out of 3490 specimens) for Microtus, Neotoma, Peromyscus, Sorex, Tamias and Thomomys were only resolved to the genus level, which we kept following previous work (Blois et al. 2010).
Abundance was reported as the minimum number of individuals (MNI), at EMC and as the number of individual specimens (NISP) at SC. At SC, MNI was additionally reported for most genus and species level observations, with exception of Thomomys sp., Sorex sp. and Scapanus latimanus. We thus used MNI at EMC and NISP at SC as inputs into metrics of functional and trait diversity.
Because the number of sampled specimens differed among layers within each site, we subsampled the layers at each site to a set number of individuals (Blois et al. 2010). At EMC, the lowest number of recorded individuals in a layer was < 10, preventing a robust calculation of community functional composition; we thus used the median number of individuals across all layers (45 individuals) as the subsampling target. All layers with fewer than 45 individuals at EMC were excluded from the analysis, leaving 14 layers, dated from 4.04 to 19.75 kybp. At SC, we subsampled to the lowest number of individuals (i.e. 130) recorded in a layer. Subsampling was repeated 500 times. Sensitivity analyses showed trends in community functional composition were robust with respect to the differences in the taxonomic resolution, the type of abundance measure and sample size (Supporting information).

Regional species pool
To produce null models (see Tests of significant change in mammalian functional diversity through time) and better define the available trait space, we expanded our assemblage to include species present in the regional species pool and their traits. We downloaded range maps of modern distributions from the IUCN Red List (IUCN 2021) for all mammals found in Spain and the United States. We defined our regional species pool as those species found within a grid cell that extended ± 1° latitude and ± 1° longitude from each site. To ensure capture of all species found in both cave deposits while still restricting the regional species pool to those species typically considered small mammals, all bats and mammals > 1100 g were removed from the dataset. We then expanded the regional species pool to include species that were present in the fossil record at these sites but have since been locally extirpated (e.g. Alborimus albipes in the SC regional species pool). As such, the regional species pool for each site reflected both modern and past potential small mammal species composition.

Mammalian functional diversity and traits
We based estimates of all functional and trait diversity metrics on a compilation of function-relevant (i.e. reflecting species' functional role in an ecosystem) traits in Wilman et al. (2014), which provides species mean values for four trait categories: body mass, diet, foraging and activity time. The diet and activity time categories included respectively ten and three axes: proportions of fruit, seeds, nectar, other plant material (e.g. leaves, described as 'plant diet' in the text), invertebrates, carrion, vertebrate ectotherms, vertebrate endotherms, fish and unknown vertebrates in species' diet (diet); diurnal, nocturnal and crepuscular (activity time).
While species can have multiple diet and activity time traits (e.g. be both diurnal and crepuscular), they are only assigned one foraging stratum: ground level, arboreal, scansorial (i.e. both ground level and arboreal foraging), aerial and marine, although the latter two were not present at either site. For all genus-level occurrences, genus-level mean trait values were used.
We used R ver. 3.6.0 (<www.r-project.org>) for all statistical analyses. Using the 'dbFD' function in the 'FD' package ver. 1.0-12 (Laliberté et al. 2014), we calculated the following functional diversity indices: functional richness (FR), functional evenness (FE), functional dispersion (FDS) and functional divergence (FDV). Functional richness describes the amount of trait space occupied by each community based on the minimum convex hull that includes all species (Villéger et al. 2008). Functional evenness describes how evenly species and their abundances are distributed in trait space such that communities with high FE are relatively equally distanced in their position in trait space, whereas communities with low FE show clustering of many functionally similar species and a few that are very functionally dissimilar (Villéger et al. 2008). Functional dispersion represents the mean of the abundance-weighted distance of species within each community from the centroid of the multidimensional trait space (Laliberte and Legendre 2010). High and low FDS values indicate that many abundant species are positioned far and close, respectively, from the centroid of the community's trait space. Finally, functional divergence measures the mean (non-abundance weighted) distance between the community's centroid and each species and then calculates the abundance-weighted distance of each species from the centroid, where species closer to the centroid are given negative distances and species with more extreme trait values are given positive distances (Villéger et al. 2008). High FDV values indicate the community has many abundant species with trait values that are far from the community's centroid, and low FDV indicates that the most abundant species are positioned close to the community's centroid.
To obtain all indices, we first conducted the principal coordinate analysis (PCoA) using the multivariate trait dissimilarity matrix based on Gower's distance for each pairwise combination of all species in the regional species pool at each site. Equal weights were given to each of the trait categories and to each axis within the trait categories (i.e. each diet and activity time axis were given a 1/10 and 1/3 weight, respectively, for a total weight of 1 for each category; body mass and foraging stratum were given a weight of 1). We used Gower's distance because it can handle quantitative, semiquantitative and qualitative variables and assign different weights to individual traits (Pavoine et al. 2009). Relative species abundance (given as MNI at EMC site and NISP at SC site) in each community was used to weight FE, FDS and FDV indices. Due to computational limitations, we based FR estimation on the first three principal coordinate axes, which preserved 64.26% and 69.85% of the total variation at EMC and SC, respectively.
We further assessed individual components of mammalian trait space -body mass, diet, substrate use and activity time -by quantifying community weighted mean trait values (Lavorel et al. 2008). All indices were calculated for each of the subsampled representations of a community, resulting in a distribution of 500 values for each index.

Temporal change in mammalian functional diversity and traits
We tested for trends in all functional and trait diversity metrics through time by regressing each metric against time using two models: one that followed a linear relationship µ i = β 0 + β 1 x i , and one that followed a quadratic relation- To select the best fitting model (quadratic or linear), we estimated the Watanabe-Akaike information criterion (WAIC), which is preferable to similar alternatives because it averages over the posterior distribution rather than using only point estimates (Gelman et al. 2014). A threshold of ∆WAIC > 7 was used to determine when model fit was decisively better (Burnham and Anderson 2004); when ∆WAIC < 7, the simpler model was chosen as the best fitting model. For each metric, we fitted 500 linear and 500 quadratic models -one for each subsampled value -resulting in a distribution of 500 intercepts and slopes for each model. The linear model was preferred to the quadratic for all functional diversity metrics and community-weighted mean trait values (Supporting information). FR and FDS are bounded by zero but have no upper bound and thus followed a gamma distribution gamma(α,β); FE, FDV and all community weighted mean trait values (with exception of body mass) are all bounded by (0,1) and thus followed a beta distribution beta(α,β). Community weighted body mass was log-transformed and thus followed a Gaussian distribution N(µ,σ 2 ). Each model was fitted in a Bayesian framework using integrated nested Laplace approximation (INLA; Rue et al. 2009) through the R-INLA package ver. 20.03.17. Priors followed a Gaussian distribution N m s , 1 2 ae è ç ö ø ÷ , with mean µ = 0 and precision 1/σ 2 = 0 for the intercept, and µ = 0 and precision 1/σ 2 = 0.001 for the slope.

Tests of significant change in mammalian functional diversity through time
For each index of functional diversity and community weighted trait values, we considered the trend through time to be significantly different from zero if the 95% credible intervals of at least 251 (50% + 1) slope estimates did not overlap zero, and had the same direction (i.e. at least 251 slopes were either positive or negative). Additionally, we statistically controlled for the association of functional diversity indices (FR, FDS, FE and FDV) with species richness by exploring whether the observed temporal trends in functional diversity indices deviated from the null expectation. Note that community-weighted mean trait values typically do not show close association with species richness, and thus do not require a statistical control (Swenson 2014). For each of the 500 subsampled communities at each site, we employed a randomization procedure wherein, for every community except the oldest, we randomized (100 times) the abundances across species in the regional species pool. For each site, we kept the relative abundance of all species in the oldest community constant in order to establish a baseline against which community change can be compared. This procedure resulted in 50 000 null expectations for each of the diversity indices (FR, FDS, FE and FDV) for each site. We tested for deviation from the null expectation in two ways. First, we calculated standardized effect size (SES) as: where 'obs' is the median observed slope value and 'exp' is the distribution of median slope values of the corresponding null communities (i.e. the expectation; (Swenson 2014)). SES values inform on a relative deviation from the null expectation; negative and positive values of SES indicate that an observed value is lower and higher, respectively, than the average expected value. We calculated SES values for all 500 observed trends (per metric and site) and report their mean values here. We further calculated quantile scores for the observed values of all indices. The observed values falling outside the 2.5% and 97.5% quantiles of the null distribution were considered statistically significantly lower and higher, respectively, than expected given change in species richness. While the quantile scores and the associated p-values remove some of the potential biases resulting from the shape of the null distribution, they also remove information regarding the effect size itself (Swenson 2014). We thus report both the SES and the quantile scores with associated p-values.

Environmental correlates of change in mammalian functional diversity through time
To understand environmental correlations with changes in mammalian functional diversity, we used climate, bioclimate and vegetation reconstructions for the last 120 000 years at 0.5° spatial resolution (Beyer et al. 2020). We extracted mean annual temperature, annual precipitation, leaf area index (LAI), net primary productivity (NPP) and biome reconstructions for EMC and SC at 1000-year intervals for the last 20 000 years (Supporting information). Because all environmental variables, with exception of leaf area index, showed significant (p < 0.05) correlations with one another (Supporting information), we transformed them into z-scores and summarized using a principal components analysis (PCA).
We tested for the relationships between each functional diversity index and community-weighted mean trait value and environmental variables using two models: one that included additive effects of PC1 and PC2 (i.e. µ i = β 0 + β 1 PC1 i + β 2 PC2 i ) and one that included an interaction effect of PC1 and PC2 (i.e. µ i = β 0 + β 1 PC1 i × PC2 i ).
For each metric, we fitted 500 models with additive and 500 models with interaction effects -one for each subsampled value -and used WAIC to compare their performance. The model that captured an additive effect of PC1 and PC2 was superior for all metrics. Using the additive model, we then tested whether the observed relationships between functional diversity indices (FR, FDS, FE and FDV) and PC1 and PC2 deviated from the null expectation given species richness, following the procedure outlined in Tests of significant change in mammalian functional diversity through time.

Temporal change in mammalian functional diversity
El Mirón Cave (EMC) experienced increases through time in both functional richness (FR; median slope = 1.85 [95% CI = 1.14-2.38]) and functional dispersion (FDS; median slope = 0.86 [0.78-0.95]) (Fig. 1). These increases were greater than change in species richness at that site would suggest (mean SES scores for FR and FDS of 4.78 . The quantile scores and associated estimated p-values did not support a significant deviation from null expectation for any of the metrics at the SC site (Supporting information). We further tested whether the observed trends are robust with regards to sample size, taxonomic resolution and relative abundance measure, and found no significant effects of these factors on the resultant functional diversity trends (Supporting information).

Ordination of environmental variables
The first two components of the PCA on environmental variables accounted for 97.6% and 97.4% of the variation at EMC and SC, respectively (Supporting information). At EMC, mean annual temperature (0.605) and NPP (0.5) loaded strongly positively, while annual precipitation (−0.594) loaded strongly negatively on PC1. LAI loaded strongly positively (0.839) on PC2. Likewise, at SC, mean annual temperature (0.589) and NPP (0.581) loaded strongly positively and annual precipitation (−0.515) loaded strongly negatively on PC1, while LAI loaded strongly negatively (−0.836) on PC2.
The SES values as well as the quantile scores and associated estimated p-values supported a significant Observed values (points) and their 75% quantiles (whiskers; based on 500 subsampled communities) are shown for dietary, activity pattern, foraging stratum and body size traits; fitted trends are shown with 95% credible intervals. Asterisks indicate that the majority of fitted trends are significantly different from zero. All dependent variables (with exception of body mass) were back-transformed to their original scales before being plotted.  . Observed values (points) and their 75% quantiles (whiskers; based on 500 subsampled communities) are shown for functional richness (FR) and functional dispersion (FDS); fitted relationships are shown with 95% credible intervals. All relationships shown are statistically significant (i.e. the majority of slopes are significantly different from zero). Mean annual temperature and net primary productivity (NPP) loaded strongly positively and annual precipitation loaded strongly negatively on PC 1 at both locations. All dependent variables were back-transformed to their original scales before being plotted. PC1 ( Fig. 5; Supporting information). At SC, invertebrate diet (median slope = 0.14 [0.12-0.17]), ground-level foraging (median slope = 0.31 [0.14-0.47]) and nocturnal activity

Discussion
We identify significant and opposing temporal trends across the late Quaternary in functional and trait diversity of European and North American small mammal communities, with differing consequences for ecosystem functioning. Figure 5. Relationships between change in small mammal community weighted mean trait values and environmental characteristics (given by principal components 1 and 2; PC 1 and PC 2) at El Mirón Cave (left) and Samwell Cave (right). Observed values (points) and their 75% quantiles (whiskers; based on 500 subsampled communities) are shown for dietary, activity pattern, foraging stratum and body size traits; fitted trends are shown with 95% credible intervals. Only statistically significant (i.e. the majority of slopes are significantly different from zero) relationships are shown. The only statistically significant relationship with PC 2 was found for body mass at Samwell Cave and is shown in grey. All dependent variables (with exception of body mass) were back-transformed to their original scales before being plotted.
In Europe, small mammal communities occupied increasingly greater functional trait space toward the late Holocene, which could not have been attributed solely to increases in species richness or an addition of one or only a few species with relatively extreme trait values to the communities. Small mammal species expanding their distributions into the region brought disproportionately unique functional characteristics that were absent in the late Pleistocene or early Holocene communities, resulting in increasing functional richness and dispersion. Specifically, through time, small mammal communities at El Mirón Cave became on average more arboreal and consumed proportionally more fruits and seeds but less of other types of plant material (e.g. leaves). This change in dietary preference suggests that the role of small mammals as vegetation engineers shifted from mostly vegetation suppression toward a combination of vegetation suppression and seed dispersal (Weltzin et al. 1997, Bricker et al. 2010. Simultaneously, the communities became less diurnal, suggesting a reduction in the amount of time small mammals spent performing other functions. In contrast, the functional trait space occupied by small mammal communities in North America became more restricted through time, and this change was closely associated with the corresponding decline in species richness (Blois et al. 2010) rather than disproportionate loss of functionally unique species. Increases in ground level foraging and insectivory paired with decreases in arboreal herbivory and declines in community mean body mass suggest that the role of small mammals in North American food webs shifted to that of secondary consumers, potentially altering the energy flow in the ecosystem (Terry and Rowe 2015).
Climatic and vegetation factors showed strong, though opposing, influences on multivariate functional diversity of small mammals at El Mirón and Samwell Caves. Increasing temperatures and slightly declining precipitation, as well as increasing net primary productivity (which, in itself, is strongly correlated with temperature) around El Mirón Cave were associated with increases in small mammal functional richness and dispersion at that location. In contrast, increasing temperature, precipitation and NPP in the region of Samwell Cave led to decreases in functional richness and dispersion, though these declines in functional diversity did not deviate from the expectation given decreasing species richness.
Our findings are consistent with climatic, vegetation and species distributions reconstructions for each region for the late Quaternary. Small mammal species of mixed Mediterranean and temperate climates, present in the Iberian peninsula during the late Pleistocene, experienced dramatic range shifts out of the region to their present habitats when temperatures warmed (e.g. Microtus oeconomus), while others immigrated into Cantabria (e.g. Cleithronomys glareolus) (Bañuls-Cardona et al. 2014). Additionally, the most recent biome reconstructions (Beyer et al. 2020) support the transition from steppe tundra to conifer and, ultimately, mixed forest in the region of El Mirón Cave around 17-15 kybp and 13 kybp, respectively (though other biome reconstructions place these transitions earlier -i.e. in late Pleistocene and early Holocene; (Carrión et al. 2010)). Such forest expansion potentially created corridors for species previously restricted to the Mediterranean climates of southern Iberia to migrate north as temperature warmed, increasing functional diversity of the region.
At Samwell Cave, the declines in functional diversity of small mammals we uncovered are likely a consequence of contractions in geographic range sizes of these species as a result of climate change (Hadly et al. 2009). Seemingly ubiquitous declines in small mammal taxonomic richness across multiple sites on the North American continent have been extensively reported by others (Graham 1976, Grayson 2000, Carrasco et al. 2009, Blois et al. 2010) and attributed to changes in climatic conditions (Blois et al. 2010, Terry et al. 2011. Despite these opposing trends in functional diversity of small mammals and their contrasting associations with environmental variables, similar environmental factors emerge as potential drivers of community change on both continents through their impacts on individual traits. Increasing temperatures and net primary productivity were associated with increasingly nocturnal activity and decreasing body mass and plant matter diet consistently across both continents. Likewise, precipitation showed consistent negative relationship with nocturnality and positive relationship with body mass and plant matter diet, despite the fact that precipitation itself showed different temporal trends at EMC and SC. Positive effects of temperature on nocturnality are consistent with the notion that small mammals respond rapidly to differences in environmental conditions and use nocturnal activity to reduce exposure to high daytime temperatures (Bennie et al. 2014). Likewise, the relationship between body mass and climate is consistent with Bergmann's rule (Bergmann 1847), wherein closely related species of homeothermic animals within a broadly distributed taxonomic clade tend to display a larger size in colder climates than in warmer regions, and corroborates findings from other studies on birds (Meiri andDayan 2003, Salewski andWatt 2017) and mammals (Meiri andDayan 2003, Lawing et al. 2017). Furthermore, higher NPP is often associated with higher taxonomic and functional diversity (Gillman et al. 2015, Williams et al. 2017) and structural complexity (Gough et al. 2019) of vegetation, likely supporting a wider array of dietary characteristics among members of small mammal communities. Both El Mirón and Samwell Caves underwent increases in NPP through the late Quaternary, with El Mirón seeing increasing forest cover and Samwell Cave transitioning from temperate grassland to open sclerophyll woodland (Beyer et al. 2020). Such increasing complexity of vegetative cover at both sites likely created more suitable habitats for arboreal granivores and frugivores (El Mirón) and insectivores (Samwell) at the expense of diets relying more on plant matter at both locations and, likely, across the larger regions of northern Spain and California.
Factors other than climate and vegetation change could have contributed to trends in small mammal community composition at El Mirón and Samwell Caves. For example, megafaunal extinction in the late Quaternary is correlated with extinctions or extirpations of small mammals (Barnosky et al. 2011a) and modern exclusion studies suggest that large mammals shape small mammal communities (Keesing 1998, Young et al. 2015. Megafauna extinctions in Europe, however, occurred in at least two pulses, during the LGM and at the end-Pleistocene (11.7 kybp) (Koch andBarnosky 2006, Stuart andLister 2007), whereas in North America large-bodied mammals went extinct in the span of ~5000 years at the end-Pleistocene (Koch and Barnosky 2006). That means that large-bodied mammals were largely extinct on both continents by the end-Pleistocene and thus played an unlikely role in the trends in small mammal diversity.
We show that small mammal communities underwent significant changes in their trait space that likely led to shifts in food web dynamics and changes in the role of small mammals as vegetation engineers in both North America and Europe. Our results point to the importance of the interaction between climate change and vegetation dynamics for community functional composition. Our findings demonstrate the utility of a trait-based approach for identifying the environmental correlates of functional diversity change through time as well as gaining valuable insight into the consequences of such change to ecosystem functioning. By building a more in-depth understanding of the ecosystem level consequences of biodiversity change across the late Quaternary, we hope to better use the lessons of the past to predict how human actions may affect ecosystems in the decades and centuries to come.

Supporting information
The supporting information associated with this article is available from the online version.