Introduction

There are several non-mutually exclusive hypotheses that attempt to explain the evolution of brain size variation within mammals. Some of these relate to selection pressures that might be imposed on the individual (Milton 1981; DeCasien et al. 2017; Rosati 2017) or social level (Dunbar 1998). Others emphasise the metabolic and developmental constraints on brain size (Aiello and Wheeler 1995; Isler and van Schaik 2009). Lastly, some emphasise the effects of intelligence and culture as capable of mitigating against both selection and constraints (Sol 2009; van Woerden et al. 2012; Holekamp and Benson-Amram 2017; Muthukrishna et al. 2018). The cognitive buffer hypothesis (CBH) is perhaps most general: it posits that larger brains have an adaptive function of buffering individuals against environmental challenges (variation in rainfall and temperature, also known as seasonality) and might facilitate mitigating selection pressures and constraints by construction of novel behavioural responses (Sol 2009). These overarching hypotheses are further contextualised through a wide range of ecological variables that have been associated with mammalian brain size, including locomotor adaptation (McNab and Eisenberg 1989; Bertrand et al. 2021), diet (Mace et al. 1981), home range, and activity pattern (McNab and Eisenberg 1989; Powell et al. 2017).

The study of brain size variation can be impacted by high levels of ‘noise’ when all mammals are studied as one taxon (Powell et al. 2017; Tsuboi et al. 2018; Smaers et al. 2021; Todorov et al. 2021). Effects of certain ‘weaker’ predictors may then be masked by more pronounced ones, or effects on smaller clades may disappear when analysed in the context of all mammals. This is particularly true when inquiry focuses on whole brain size, instead of specific brain partitions (Barton and Harvey 2000). Therefore, testing different hypotheses might be best done by focusing on smaller clades and on different brain components. Lagomorpha (rabbits, hares, and pikas) are particularly suitable in this respect. This mammalian order is part of the superorder Euarchontoglires, which also contains primates (along with rodents, tree shrews, and colugos). Lagomorphs appear in the fossil record around 52 million years ago in present day China and Mongolia (Lopatin and Averianov 2008; Wang et al. 2010). The majority of lagomorphs are rabbits (42 species) and hares (33 species), together grouped in the family Leporidae, but they also include the more plesiomorphic pikas (34 species) (Murphy et al. 2001). However, leporids have a number of shared and highly uniform traits that may reduce the noise in modelling evolution of brain size.

The uniform leporid locomotor mode emphasizes adaptations for increased half-bounding cursoriality (Kraatz et al. 2021), with which some species can reach high maximal speeds (up to 72 km/h in Lepus europaeus and Lepus alleni; Garland 1983). This specialised locomotion appears to be associated with a highly unusual case of cranial kinesis, where the alisphenoid and squamosal bones – which also form part of the braincase - are separated by the piriform fenestra and can move relative to each other (Wible 2007). Cranial kinesis has been proposed to be related to the demands of locomotor behaviour, dissipating kinetic energy between the loading cycles of cursorial locomotion (Bramble 1989). In addition, all extant leporids are exclusively herbivorous and minimally sexually dimorphic, and their life history is largely adapted to high predation which is related to both their large litter sizes and abundance cycles (Krebs et al. 2001). This uniformity is also reflected in what appears to be a highly conserved anatomy of the brain. The first virtual endocast of a fossil lagomorph, Megalagus turgidus, described by López-Torres et al. (2020), despite indicating larger olfactory bulbs (OBs) compared to extant leporids, exhibits only a lateral sulcus on the neocortex on an otherwise lissencephalic cerebrum. This means that brain anatomy in the lineage is consistent throughout evolution and has been stable since the emergence of stem lagomorphs.

The striking commonalities among leporids (locomotion, cranial morphology, diet, low sexual dimorphism, and similar life history) have been associated with variation in brain size in other lineages in the past (Kraatz et al. 2015; DeCasien et al. 2017; Muchlinski et al. 2018; Weisbecker et al. 2021), making them a relevant example for assessing brain size with these factors largely controlled. At the same time, leporids inhabit diverse habitats (Kraatz et al. 2021), which is known to be associated with brain size variation (Gonda et al. 2013; Corfield et al. 2015; Sayol et al. 2018; De Meester et al. 2019). In addition, leporid reproduction is vulnerable to environmental variability (Tablado et al. 2012), which can result in high (50%) mortality within litters (Rödel et al. 2009). This makes leporids also an ideal clade to test hypotheses of a relationship between environmental variation and brain size, such as the cognitive buffer/seasonality hypothesis.

Within lagomorphs, the only species that has been studied extensively in terms of brain structure and function is the domesticated European rabbit (Oryctolagus cuniculus). This species has a long history of co-existing with humans as a food source; leporids were domesticated by humans in the 5th century in Europe by French monks (Clutton-Brock 1989), and they used to be a common model organism. However, little is known about the evolution of the lagomorph brain (Kraatz et al. 2021). Breeds of domesticated rabbit have been used in studies of brain developmental timing (Lim et al. 2015), development of the visual cortex and neuronal morphology (Murphy and Magness 1984), brain anatomy (Schneider et al. 2018), and the effects of domestication on the brain (Brusini et al. 2018). Despite the focus on the domesticated rabbit, larger scale studies on brain evolution in lagomorphs have not been performed to date. The first virtual endocast of a fossil lagomorph, Megalagus turgidus, was described by López-Torres et al. (2020). It has larger olfactory bulbs relative to its endocranial volume compared to extant leporids. The frontal lobes of Megalagus are wider than those observed in modern lagomorphs and not as expanded rostrally. Similar to extant leporids, the endocast of Megalagus exhibits only a lateral sulcus on the neocortex on an otherwise lissencephalic cerebrum, meaning brain anatomy in the lineage is consistent throughout evolution and has been stable since the emergence of stem lagomorphs around 42 Ma.

The radiation of leporids is also excellent for assessing if and how variation in the size of the olfactory bulb (OB) reflects environmental variation. One of the primary neurological systems processing chemical molecules in mammals is the OB, comprised of the main and accessory olfactory systems (Swaney and Keverne 2009; Bertrand et al. 2022). The OB is one of the earliest structures to evolve in complex central nervous systems (Eisthen 1997). It is heavily involved during sexual and social behaviours, as well as in spatial orientation (Ihara et al. 2013) and it has a different neuronal scaling relationship from the rest of the brain (Finlay et al. 2001; Ribeiro et al. 2014). Moreover, it is extremely plastic throughout ontogeny (Wu et al. 2020) and highly variable in size and structure across vertebrates (Finlay et al. 2001). Previous research has shown that the OB has undergone both reduction and expansion in different mammalian lineages, independent of the rest of the brain (Meisami and Bhatnagar 1998; Finlay et al. 2001). Additionally, the evolution of the OB is related to the evolution of thermoregulatory function in some mammals (Martinez et al. 2020). Moreover, olfactory behaviour, mediated by the neural circuit devoted to olfaction, has been shown to be important for processing environmental spatial and social information in mammals and birds (Benhamou 1989; Corfield et al. 2015).

Leporids are an ideal radiation to study variation in the olfactory bulb because they have a highly developed and conserved double path vomeronasal structure while its cellular structure is highly variable among species and between sexes (Villamayor et al. 2020). Previous studies have shown that in mammals, rabbits represent one of the best models for studying chemocommunication (González-Mariscal et al. 2016) and they remain the only mammal species in which a mammary pheromone (2-methylbut-2-enal) has been comprehensively characterised (Schaal et al. 2003). This pheromone is released by lactating females to awaken rabbit neonates and initiate the nipple-sucking reflex. Pheromones are a key means of communication in rabbits and are important in submissive and dominant behaviours (Melo and González-Mariscal 2010). Due to rabbits often being prey animals, olfaction is an essential sense for the detection of danger, predators, and mates (Apfelbach et al. 2005). This heavy reliance on olfaction in different aspects of leporid’s natural history makes the OB an excellent target for brain evolution studies in the family.

In this study, we use Bayesian phylogenetic comparative methods in combination with multiple imputation technique to test five models of brain and olfactory bulb evolution using a sample of 18 species of leporids derived from 60 individual virtual endocasts from (Kraatz et al. 2015). Our first hypothesis is related to the Cognitive buffer hypothesis, where we expect seasonality (as measured by seasonal variation in precipitation and temperature) to have a negative relationship with both OB and rest of the brain (ROB) in leporids, because of the challenges imposed by unpredictable food sources, sheltering opportunities, and predator behaviour. Alternatively, if seasonality does not affect brain size in this family, this might be due to alternative behavioural strategies like burrowing, expansion of home or geographic range. Another set of hypotheses we test are related to maternal investment, where we expect litter size and gestation length to be negatively related to brain size (Isler and van Schaik 2012; Weisbecker et al. 2015; Todorov et al. 2021), and dietary/ecological pressures as measured by diet niche breadth (Milton 1981). Previous studies have provided conflicting findings regarding the relationship between brain size variation and diet in other mammals (DeCasien et al. 2017; Todorov et al. 2019) so given that all leporids have comparable diet, we are controlling for the confounding effects of this variable.

We focus on the OBs due to their prominence and the ease of measurement from brain endocasts, in addition to their important functional involvement in odour communication and spatial orientation (Pager 1986; Gittleman 1991; Ihara et al. 2013; Yi et al. 2021). Moreover, it has been shown that the OBs in mammals do not follow the same allometric scaling relationship to body size as the ROB (Finlay et al. 2001) thus making it a good candidate for brain evolution studies.

We expect that the Olfactory Bulb and Rest of Brain in leporids will be under different selective pressures (López-Mascaraque and de Castro 2002; Treloar et al. 2010). Specifically, we expect the ROB to be more constrained by developmental factors, while the OB will be more sensitive to selection pressures, namely seasonality. Seasonality (especially in precipitation) has been shown to affect the perception of chemical cues, whereas environments with higher precipitation would not be conducive for olfactory orientation due to the reduced transmission of volatile chemical cues (Bowen and Janzen 2005; Czaja et al. 2018). On the other hand, it might be that the effects of seasonality can be buffered through behavioural and cognitive adaptations such as burrowing or activity cycle adjustment (diurnal vs. nocturnal) instead of through changes in the OB and ROB size and architectonic.

Materials and Methods

The virtual reconstructions of endocasts we used are available on MorphoSource (https://www.morphosource.org/projects/000417960 for the endocast sample and https://www.morphosource.org/projects/000415700 for the cranial sample). All code is available online on https://github.com/orlinst/Leporid-brain-evo.

Study Animals

Morphometric and ecological data were collected from 60 leporid skulls spanning 18 species within all extant genera of Leporidae. Most specimens were from the American Museum of Natural History (AMNH) or The Los Angeles County Museum of Natural History (LACM). For further detail on the origin of the skull specimens, refer to (Kraatz and Sherratt 2016) and the MorphoSource dataset. The sample reflects the phylogenetic coverage in Kraatz and Sherratt (2015), based off of the Matthee et al. (2004) leporid phylogeny. All specimens used were adults, characterised by fully fused occipital sutures (Hoffmeister and Zimmerman 1967).

Ecological Data

Ten ecological and life history variables relevant to our hypotheses were collated from existing literature including: activity cycle (diurnal, nocturnal, crepuscular), locomotor mode (generalised, saltatorial, cursorial), diet breadth (three categories as defined by the number of host plants used, categorised by PanTHERIA (Jones et al. 2009) and used as a continuous variable in our analyses), gestation length, home range, geographic range, litter size, burrowing behaviour and seasonality in precipitation and temperature (Kraatz et al. 2015; Myers et al. 2022). PanTHERIA was used to source the following variables: adult body mass (grams); diet niche breadth; gestation length (days); home range (km2); litter size, and geographic range (total extent of species range with a global equal-area projection). All continuous variables were natural log transformed. Average temperature and precipitation seasonality was compiled for each leporid species within its range. For species with contiguous geographic distributions, at least 10 years of data for monthly temperature and precipitation were sampled from KNMI Climate Explorer database (climexp.knmi.nl) using every available weather station within minimum and maximum latitudes and longitudes within each species geographic range. For species with discontiguous ranges or low/no weather stations in the KNMI Climate Explorer database, monthly temperature and precipitation were sampled randomly from weather stations within the species geographic range from Weatherbase (www.weatherbase.com). The average annual seasonality for combined years and all sampled weather stations for both temperature and precipitation within a given species range was defined as the coefficient of variation:

$$S=\frac{SD}{M}$$

Where within the geographic range of a given species, S is the annual seasonality of temperature or precipitation, SD is the standard deviation of mean monthly temperature (oC) or precipitation (mm), and M is the mean monthly temperature (oC) or precipitation (mm).

Temperature values were first converted to absolute temperature by adding 273.15 and, after using the above formula, the results were multiplied by 100.

3D Reconstruction of Endocasts

Virtual endocasts from skulls of the 60 leporid specimens were prepared from micro-CT scans obtained by Kraatz and Sherratt 2016; (voxel size 125-160 μm)  using Mimics version 18 (Materialise NV, Leuven, Belgium, 1992–2020). Brain endocasts were prepared through “flood-filling” the cranial cavities of scanned specimens in Mimics and 3Matics by NH.

Endocranial Volume Measurements

Measuring brain partitions from endocasts can be challenging and, given the quality of our sample, we were able to measure the size of the OB along the whole endocranial cavity (but see Carlisle et al. 2017). After measuring total brain volume, the endocasts were divided into OB volume and ROB (total minus olfactory bulb volume). To compartmentalise the brain into OB and ROB, brain endocasts were virtually segmented using the software Mimics. To delimit the posterior end of the brain consistently from the spine, the brain stem posterior to the cerebellum (and the opening of the foramen magnum) was not reconstructed. The endocast thus ended along a flat plane from the dorso-posterior cerebellum to the floor of the basicranium. This was done to standardise brain stem size between endocasts, which may otherwise confound volumetric analyses. OBs were digitally separated by the placement of a coronal surface plane immediately posterior to the proximal point of the olfactory bulb (See Online Resource 1, Section 7, Figs. 1 and 2). Partition volume for the OB and ROB endocast were calculated using Mimics. All endocranial volumes were independently measured by two researchers (CH and MM), and the measurements were identical with r = 0.998. Endocasts of all specimens used in this study are available at MorphoSource.

Fig. 1
figure 1

Phylogenetic tree of the 18 species used in the study, based upon Matthee et al. (2004) with heat map displaying all min-max normalised mean species values of continuous variables. The lightest cells represent the lowest mean values (0) per species per variable and the darkest cells represent the highest values (1). White cells represent missing data. Blue circles indicate burrowing species, and yellow squares indicate non-burrowing species. Explanation of variables: Bmass, body mass (kg); Brain, total endocranial volume (mm3); GesLen, gestation length (days); LitterS, letter size (number of offspring); Range, home range (km2); SeP, seasonality in precipitation; SeT, seasonality in temperature

Fig. 2
figure 2

ANCOVA of total brain volume plotted against body weight of all specimens in the sample (n = 60). Different species are indicated in different colours (listed in the legend in increasing order related to body size and the number of individuals per species in brackets). The red line and polygon indicate the regression line of burrowing species, and the green line and polygon represent non-burrowing species. The Xs mark the mean value per species used in the regression analyses

Statistical Analyses

All analyses were conducted in the R statistical environment v4.1.1 (R Core Team 2021) using the following packages for the analysis: MCMCglmm (Hadfield 2010), MulTree (Guillerme and Healy 2014), ape (Paradis and Schliep 2019), phytools (Revell 2012), mice (Buuren and Groothuis-Oudshoorn 2011), phylomice (Drhlic P, Blomberg SP, unpubl.; available from SPB), naniar (Tierney et al. 2019), ggplot2 (Wickham 2016), bayesplot (Gabry et al. 2019), hdrcde (Hyndman et al. 2021) and performance (Lüdecke et al. 2021) for plotting.

Phylogeny

We used a tree derived from Matthee et al. (2004) that includes the 18 species studied here (Fig. 1). The tree had two polytomies that needed to be resolved due to the requirements of some of the analyses. We did so by adding 0.01% of the median branch length. Subsequently, the tree was ultrametricised using the function ‘extension’, from the R package phytools.

Imputation

Testing hypotheses related to brain evolution within small clades inadvertently suffers from issues related to sample size, constraining the scope of such inquiries. Issues related to limited sample sizes include low statistical power, faulty estimation of various model parameters, and limited scope of the analytical conclusions. Because of the scarcity of data, both of anatomical traits and traits related to life-history and ecology, many studies omit whole cases due to one or a few missing traits, as most statistical comparative methods cannot deal with missing data. One way to solve this pervasive issue is through multiple imputation (Rubin 1987). In this study, we have used phylogenetic multiple imputation to conserve the complete sample (n = 18) for all analyses.

Our dataset contained 3.3% missing data (with four/22% of species missing Home Range or Diet Breadth values, and one/5.5% of species missing a Gestation Length value). No phylogenetic signal in the missingness pattern was detected, as shown by our analysis using the R package sensiPhy with D statistic (Fritz and Purvis 2010). Analysis of the pattern of missingness are included in the data available at https://github.com/orlinst/Leporid-brain-evo (see also Online Resource 1) and have been conducted using the R package naniar (Tierney et al. 2019).

We imputed four datasets using the phylogenetic predictive mean matching algorithm in phylomice, and all four imputed sets were used in all models including variables with missing values (Home Range, Diet Breadth, and Gestation Length). For details regarding the imputation algorithm and protocol see Rubin (1987), Little ( 1988), Barnard and Rubin (1999), White et al. (2011), and Todorov et al. (2021).

Modelling

We tested five different models as described in Table 1.

Table 1 Variables included in the models. Each model included body mass as a covariate, and olfactory bulb (OB) or rest of the brain (ROB) were used as dependent variables. Models containing variables with missing data (marked with asterisk) were run on multiple imputed datasets

All models were analysed using MCMCglmm (imputed datasets) and pGLS (complete cases with case-wise deletion included in the supplement; see Online Resource 1). A full description of the models tested is included in Online Resource 1. Each model was assessed using the performance R package, and all met the assumptions for linear regression (see Online Resource 1).

For the MCMCglmm models, we used improper, uniform, uninformative priors, which assume that all values of the parameters are equally likely. We ran the MCMCglmm for 1,000,042 iterations with a burn-in of the first 100,000 iterations and a sampling rate of 300. Each model was run on two chains and resulted in effective sample size of at least 2,000. All models converged successfully using the Gelman-Rubin criterion < 1.1 (Brooks and Gelman 1998).

The results from models based on the multiply-imputed datasets (two parallel chains of MCMC on four imputed data sets) were pooled using Rubin’s rules (Rubin 1987). The fit of all models was compared using the deviance information criterion (DIC) (Spiegelhalter et al. 2014). This is an estimator of prediction error, similar to AIC, where the estimate is based on the posterior mean. Only model parameters with substantial posterior distribution away from 0, defined as at least 95% above or below 0, were selected as being significant.

For the pGLS analyses, we used the R package ape, using Pagel’s correlation structure and omitting all cases with missing values. We conducted pGLS analyses on all models tested with MCMCglmm and on litter size and seasonality following up on a point raised in the discussion.

Results

We set out to test the effect of seasonality on leporid brain evolution and, more specifically, on the variation of OB and ROB volumes. We did not find any support for such an effect in either case (see Table 2). We tested four other models related to hypotheses explaining the evolution of brain size variation (shown in Tables 1 and 2), and we found support for two of them: leporids with larger home ranges that exhibit burrowing behaviour have larger brains than non-burrowers (bolded in Table 2) or species with smaller home ranges. Species with larger home ranges were shown to have both larger OB and ROB, while burrowing species were shown to have larger ROB and marginally significant larger OB (94% of the probability distribution above 0, with a threshold of 95%). Our pGLS analysis (Online Resource 1) included some Pagel’s lambda values (indicating the strength of phylogenetic signal in the residuals) that were negative or higher than 1, which violates the definition of lambda per se (see Discussion).

Table 2 Results from the MCMCglmm analysis on the full MI dataset. Abbreviations: ROB, Rest of brain volume; OB, Olfactory bulb volume; β, regression coefficient; SE, standard error; PD > 0, posterior distribution above 0 in percentage; DIC, deviance information criterion

Discussion

We set out to test five different hypotheses related to brain evolution in leporids: we expected seasonality (as measured by seasonal variation in precipitation and temperature) to have negative effects on both OB and ROB or, alternatively, if seasonality does not affect brain size, alternative behavioural strategies like burrowing or expansion of home or geographic range were expected to relate to an increase in brain size. Another set of hypotheses we tested were related to maternal investment, where we expected litter size and gestation length to be negatively related to brain size, so as dietary/ecological pressures as measured by diet breadth.

Unexpectedly, we found no support for any effect of seasonality on the evolution of brain size variation in leporids. On the other hand, our analysis indicated that the only correlates of brain size in this family are burrowing behaviour and home range. Seasonality in temperature or precipitation has previously been shown to have both positive (van Woerden et al. 2012) and negative (Luo et al. 2017; Heldstab et al. 2018) effects on the evolution of brain size in vertebrates, but we were unable to detect any effect in our sample. One reason for the lack of such effect in leporids might be due to the fact that burrowing behaviour compensates for climatic variation, and as such, could be a selection force responsible for an increase in brain size. Another possible mechanism for overcoming seasonality pressures due to unpredictability of food sources and predatory pressure might be through extending home ranges which is also supported by our findings. This is in line with the cognitive buffering hypothesis (Sol 2009) where individuals can overcome some environmental challenges by behavioural and cognitive adaptations. Alternatively, seasonality may be affecting other aspects of brain evolution in leporids, besides brain volume (i.e. architectonic, physiological, morphological) for which further investigation is needed.

Our analyses do not support the notion that species with larger geographic ranges are better at avoiding seasonality selective pressures related to brain size. The only variable that is positively related to brain size from our spatial model is the species’ home range. Despite leporid home ranges being relatively small (from 0.004 to 0.86 km2) larger home ranges might provide better opportunity for niche construction through burrowing, predator avoidance and other opportunities related to seasonality buffering (Bond et al. 2001).

The fact that the size of the olfactory bulb was shown to be related to the size of the home range in leporids is consistent with the finding in mice that processing of olfactory information is related to processing of spatial information (Wiedenmayer et al. 2000). It has also been previously shown that the size of the olfactory bulb is related to home range in carnivores (Gittleman 1991) and this confirms that larger olfactory bulbs are beneficial in maintaining extended home ranges and mitigation of predation and seasonality pressures. However, the exact mechanism behind this relationship in not yet clear.

We did not find support for any maternal investment effect (as measured by litter size and weaning age) limiting brain size evolution in leporids. This is unexpected because many previous studies have found that litter size correlates negatively with brain size. This might be due to the fact that our sample is relatively homogeneous in terms of these two reproductive variables, and while such an effect has been shown in larger clades (marsupials (Todorov et al. 2021) and mammals in general (Isler and van Schaik 2012)), focusing on the family level, such an effect might not be detectable in our small sample. Additionally, the limiting effect of maternal investment on brain size variation might be ubiquitous and uniform in most leporids, as in the lineage (and in other animals), larger litter sizes are maintained to buffer environmental variability, predation and infectious disease (Bond et al. 2001; Reddon et al. 2018; Vega-Trejo 2021). It has previously been shown that an increase in seasonality is related to an increase in litter size in rodents (Stewart et al. 2020). We confirmed this in our sample of leporids (pGLS Litter Size ~ Seasonality in Temperature p = 0.0336, beta = 0.66, t = 2.32, λ = 1.07, df = 16; but no relationship to Seasonality in Precipitation (p = 0.64, t = -0.48, df = 16). This might indicate that one of the evolutionary mechanisms buffer seasonality in rodents and leporids is reproductive rate, rather than cognition.

Additionally, our full activity model was also not supported (including locomotor mode and activity cycle) indicating that leporid species, being mainly terrestrial, are unable to buffer any environmental effects relating to brain size, besides using burrowing behaviour, and thus differences in activity mode are not related to brain variation.

Surprisingly, in most of our pGLS models (included in Online Resource 1) including ROB and in one including OB, the lambda estimates were unreasonable. This is a common situation in many phylogenetic comparative studies (De Meester et al. 2019) but it is often ignored. Pagel’s lambda (Pagel 1999), as a measure of phylogenetic signal in the residuals, can only take values between 0 (indicating lack of phylogenetic signal), and 1 (indicating fit with Brownian motion), so any value larger than 1 or lower than 0 is a result of a flawed restricted maximum likelihood estimation and therefore unrealistic. However, in all pGLS models including ROB, the estimated maximum likelihood lambda values always fell outside the 0–1 interval (see Online Resource 1). This is most probably due to the small sample size, which renders all our lambda estimates unreliable. This and other cases of the same issue (De Meester et al. 2019) warrant caution when using small sample sizes in phylogenetic comparative studies, as then, maximum likelihood estimations become unreliable. Unfortunately, the limitations of sample size are common in many comparative studies, and while such analyses lay a foundation for further inquiries after more extensive data collection, they should be taken with substantial caveats. Authors of such studies should pay more attention to the statistical artefacts arising from small sample size limitation and analyse them in more detail before reporting their results (Freckleton 2009). One potential solution to this problem is doing the analysis assuming Brownian motion (λ = 1) and verifying whether the outcome of the regression with the flawed estimate corresponds to the one with the fixed lambda value.

Our study has several limitations, most evidently relating to the common issue of low of specimen availability for comparative brain size studies, and associated unrealistic estimation of Pagel’s λ. Although covering 25% of leporid genera with the specimens available to us and having good within-species sampling of endocranial size, a larger sample would have allowed the use of more refined models. Given the conserved body plan and the relative uniformity of life-history within the family, we do not expect the inclusion of new specimens to drastically change the results of the current study, but it will improve the statistical power and make some marginal effects (i.e., the relationship between OB and burrowing) clearer. In addition, the extensive field observations required for behavioural data means that good coverage is difficult to obtain at a large scale for most mammals, such that some of our models required imputation. With increasing availability of virtual data (Davies et al. 2017), studies like ours will hopefully have better specimen coverage in the future. Nevertheless, our results confirm that consideration of a specific clade of mammals has potential to produce a much finer-grained and biologically realistic perspective on the drivers of brain size evolution.