Distribution modelling of the garden dormouse Eliomys quercinus (Linnaeus, 1766) with novel climate change indicators

The garden dormouse Eliomys quercinus has been declining in both abundance and range since the mid-twentieth century. The eastern edge of its range has contracted from the Ural Mountains to eastern Germany. Habitat loss and fragmentation has been the most supported theory to explain the observed decline. Climate change has been implicated in declines of other terrestrial mammals, but not investigated for E. quercinus. To better understand the factors influencing the distribution of this species and to map habitat suitability for E. quercinus across Europe, we created a Maxent species distribution model. Among the main environmental variables used for the modelling, two novel climate change indicator variables were produced to indicate the degree of climate change between the early twentieth century and the present. Areas of high suitability were mapped, and variable importance estimated using jackknife tests and variable contribution metrics. The climate change indicators outperformed many conventional variables, which could indicate that climate change is a factor behind the current distribution of E. quercinus. We also analysed the land use types where presence points of E. quercinus were located and whether they were in areas of “high nature value farmland”. Over 30% of all spatially filtered presence points corresponded to high nature value farmland areas. Our results could indicate a role for changing climate (particularly in temperature) in the range decline E. quercinus, and for high nature value farmland practices in conserving this species. Field studies and improved monitoring for this species are recommended to confirm both possible findings.


Introduction
The earth is undergoing its 6th mass extinction, with countless species of flora and fauna lost and many at risk of extinction in the near future (Ceballos et al. 2017). Many anthropogenic drivers have been identified behind these losses such as climate change (Morrison et al. 2018), habitat loss and fragmentation (Staude et al. 2018;Mazgajski et al. 2010), direct human-wildlife conflict (Buchner et al. 2018), the introduction of invasive alien species (Václavík and Meentemeyer 2012;MacDonald and Harrington 2003), agricultural intensification (Bertolino 2017) and pollution such as pesticide use (Schreinemachers and Tipraqsa 2012). These anthropogenic pressures (particularly habitat fragmentation) have also effected small mammals in many cases (Bertolini 2017;Palmeirim et al. 2018) and invasive species (Zozzoli et al. 2018), although climate change is an emerging concern (Williams et al. 2015). Defaunation (the extinction or extirpation of animal species) can also affect ecosystem function, and further impact other species.
Anthropogenic climate change is the change in climate conditions driven by human activities (Opdam and Wascher 2004) with significant effects including changes to temperature, precipitation and seasonal patterns, and higher frequency of severe climate events (including storms, droughts, Handling editor: Mauro Lucherini. and extreme winters) (Morán-ordóñez et al. 2018). Climate change has also increased the impact of other biodiversity threats with a synergistic effect. For example, it has facilitated the spread of invasive species (Hellmann et al. 2008) and damaged species populations (through chronic climatic shifts beyond environmental tolerances) or by impacting wider ecosystem functions). Climate change has also increased the number of severe climatic events (such as severe unseasonal storms) while habitat fragmentation may prevent rescue or recolonization of populations decimated by climate change (or a severe climate event) (Opdam and Wascher 2004).
Among the declining mammalian species is the garden dormouse Eliomys quercinus. E. quercinus is a small, primarily arboreal rodent of the order Gliridae that is currently classified as "Near Threatened" by the IUCN (Bertolino et al. 2008). E. quercinus was originally spread across almost all of continental Europe ranging from Iberia in Portugal at the west to the Russian Ural Mountains in the east (Filippucci 1999;Ruiz and Román 1999). From the mid-twentieth century, starting in the east of Europe, it began to experience declines and various populations have been completely extirpated. Indeed, Bertolino (2017) estimated E. quercinus has lost over 60% of its original range and has now been driven back to western Germany. Several reasons for this have been proposed, such as land use change (Bertolino 2017), competition with the brown rat (Rattus norvegicus) (Nowak 1996) and in some cases, human-wildlife conflict (Buchner et al. 2018); for example, some homeowners have been known to poison E. quercinus that live in their homes and there is a long tradition of hunting dormice species in Italy (Carpaneto and Cristaldi 1994). However, in many areas, these pressures alone seem inadequate to justify the disappearance of the species (for example, significant areas of relatively unfragmented forest still exist in eastern Europe where E. quercinus has been extirpated) and E. quercinus is also known to thrive in pine plantations (Benayas et al. 2017), which are abundant in eastern Europe. E. quercinus was last observed in north eastern Germany in 1977 (Borkenstein 1993) but is still present throughout southern and western Germany (Meinig and Buchner 2012) and its range decline has apparently slowed, perhaps as a result of improved conservation measures or local factors. It also remains unclear why E. quercinus has undergone a more extreme range contraction than other similar dormouse species. For example, Glis glis is actually abundant and increasing across Germany and is surviving as far east as Lithuania (Juškaitis et al. 2018). Nonetheless, landscape factors such as land cover or the density of roads (which could indicate habitat fragmentation or exert roadkill pressure) could have had important effects on E. quercinus' distribution.
E. quercinus provides ecosystem functions, such as seed dispersal and invertebrate predation (Bertolino 2016a) and is also a prey for several species such as sparrowhawks (Accitpiter nissus). Bertolino (2006) asserted that habitat fragmentation and agricultural intensification were the most likely causes behind the species range loss. Therefore, both land cover and land management changes could have driven changes in the species' distribution. High nature value farmland (HNVF) areas could have been beneficial to many European species (Bernués et al. 2016), including E. quercinus. Indeed, HNVF have generally been managed with less intensive use of fertilisers and pesticides than most agricultural lands and preserve features like hedgerows and meadows.
Species distribution modelling (SDM) methods, also called ecological niche modelling (ENM), are correlative modelling techniques which have become prominent in ecological study. One of the most widely used SDM methods has been Maxent (Philipps and Dudık 2008), a presencebackground-based approach which models the habitat suitability for a species based on observations of presence and environmental variables such as temperature, precipitation, and elevation (Philipps and Dudık 2008). Correlative approaches have been shown to be particularly effective when there are abundant presence records available (Elith et al. 2010). These techniques have been used to estimate potential distributions of small mammals and identify conservation priorities (Prieto-Torres and Pinilla-Buitrago 2017) as well as to guide field surveys when resources are limited or data are unclear (West et al. 2016). Moreover, modelling techniques could be adapted to predict potential suitable habitat in currently poorly surveyed areas. Presence-background only methods have some limitations when compared to presence-absence methods-in particular, these methods have struggled with unknown or variable prevalence (the proportion of suitable sites actually occupied by the target species) and the effects of sampling bias (differing levels of sampling effort or detection probability across the landscape) (Fithian et al. 2015;Stolar and Nielsen 2015). Correcting or adjusting for these issues has thus proven essential to ensure appropriate, reliable outputs although the problems cannot be completely cured (Guillera-Arroita et al. 2015). The primary alternative, presence-absence models, can avoid these issues but they usually require dedicated sampling studies to gain true absence data (thus rendering them impractical for a wide-ranging species such as E. quercinus).
The threat of future climate change could overshadow the effect of climate change which has already occurred in recent decades; across multiple taxa and geographic scales (Parmesan 2006;Stewart et al. 2015). Uribe-Riveira et al. (2017), did consider past climate and historical presence records when building their models, by calibrating their models using climate data for the period 1950-1970 to model the distribution of Darwin's frog Rhinoderma darwinii in South America and then projecting the model to the present, although the rate of climate change within their study region was not considered as a variable. Although recently occurred climate changes have not commonly been considered in distribution models, combining conventional variables with historical climate change indicators (e.g., indicators of a climate change that occurred from the beginning of the twentieth century) could contribute to better understand species distribution. Given the sensitivity of dormice to thermal stress (Santoro et al. 2017) and the known effects of winter climate change on many species (Williams et al. 2015;Inouye et al. 2000), then the rate of climate change on extreme minimum and maximum temperatures could be relevant factors behind their remaining distribution.
The objectives of this study were to identify the main environmental variables (including newly created climate change indicators) influencing the current distribution of E. quercinus and map the areas of relatively high habitat suitability for this small rodent.

Origin of data
Data on E. quercinus presence were gathered from the Global Biodiversity Information Facility (GBIF 2020), Vertnet (Vertnet 2020) and the national species presence databases of Italy (Network Nazionale Biodiversita 2018), Finland (Finnish Biodiversity Information Facility 2020), and Switzerland (Swiss Center for Cartography of Fauna 2018). We filtered the presence point data to include only those points observed between 2000 and the present, because E. quercinus may become extinct in areas where it was recorded before year 2000 (for example, presence points exist for Estonia from the 1980's, but the species has since been extirpated from that nation).
Where possible, the most temporally representative reliable data were used for each variable. Data on general land use for observation proportions were taken from the Copernicus Project, using the CORINE 2012 dataset (Copernicus Project 2020). The refined 2006 data for high nature value farmland (HNVF) were obtained from the European Environment Agency (Paracchini et al. 2008). In addition, we used the 19 standard bioclimatic environmental variables (using Worldclim version 2) with 30 s resolution (Fick and Hijmans 2017), GMTED2010 elevation (DEM) and slope data from Amatulli et al. (2018), land use across the whole training and prediction area for year 2015 from the European Space Agency Climate Change institute land cover data (ESA 2020) (the Corine data, though available at higher resolutions, only covers the western European nations within our study region, and does not cover the Ukrainian, Belarussian or Russian parts of the species former range) and finally data on the percentage cover of open water, broadleaved forest and coniferous forest (Tuanmu and Jetz 2014). All these data had a native resolution of 30 arc seconds (approximately 1 km 2 ), which was retained (both to prevent distorting the data, and because E. quercinus home ranges are smaller than 1 km 2 (7.5 ha; Bertolino 2016).

Synthesis of climate change indicators
We produced two novel variables to indicate recent historical climate change, using the CRU TS v4.03 gridded historical climate data (0.5° resolution) produced by the University of East Anglia Climate Research Unit (Harris et al. 2020). Using Arcmap v10.6, the mean maximum and mean minimum temperature data were extracted for the periods 1901-1920 and 2000-2018. The difference between these two periods was then calculated using Raster calculator in Arcmap v10.6 to produce indicators of how maximum and minimum temperatures have changed (hereafter, we refer to these variables as MinTempChange and MaxTempChange) across our study area since the beginning of the twentieth century. Finally, the data were resampled to match the 30 arc second resolution of the other variables.

Bias correction
The full database of presence points was compiled in Microsoft Excel and filtered to remove records without any geographic coordinates and records from outside the temporal focus (pre-2000). The remaining data (8214 records) were then imported into Arcmap v10.6 and spatially rarefied by 3 km using SDMToolbox (Brown et al. 2017) to reduce overfitting to sampling bias, decrease spatial autocorrelation and to remove duplicate points (Boria et al. 2014). The 3 km rarefaction was chosen as a very conservative figure, because our variables have an approximate 1 km 2 resolution and E. quercinus has a small home range. The final total number of spatial records post-rarefaction and temporal filtering for the model was 3517, the full list of points used can be found in supplementary material (Online Resource 1).
The use of bias correction has been shown to significantly improve the reliability of Maxent models (Fithian et al. 2015;Kramer-Schadt et al. 2013). In the current analysis, we considered the existing presence data to be biased due to differing sampling intensities within the different countries across E. quercinus' range (for example, Spain had many records whilst neighbouring Portugal only had two); thus, we identified target group density background (TGB)-based bias correction (Warren et al. 2014) to be the most appropriate method. To perform TGB correction, we downloaded records of occurrence for small mammal species with similar sampling methodologies (in our case, Rodentia and Soricidae) from GBIF for the entire training extent, for the period 2000-2018 and produced a heatmap to indicate sampling effort across the training area 1 3 using SDMToolbox (Brown et al. 2017). This heatmap indicates which areas have high or low sampling bias for small mammals-Maxent uses this to alter the spatial density of background points when training and evaluating the models.

Correlation of variables and stepwise refinement procedure
All predictor variables were imported to R v 3.1.5 (R core team 2019) and clipped to the study extent using the mask command in raster package. A variance inflation factor (VIF) test (we removed variables with VIF > 10) and pairwise plots were performed to determine which variables, if any, were significantly colinear: if two variables were correlated with each other, then one would be removed for sake of improving data validity (Elith et al. 2010).
After removal of correlated variables, the 10 remaining variables were: Bio6 (minimum temperature of coldest month), Bio14 (mean annual precipitation), Bio15 (annual precipitation range), land use, density of roads, slope, percentage of broadleaved forest cover, percentage of coniferous forest cover, percentage of open water and DEM (elevation).
Variable set was refined using a reiterative stepwise removal procedure similar to Zeng et al. (2016). Prototype runs were produced using our main model methodology (see "Methods" section e.), but with all variables. Variables with low contribution (< 10%) were then removed in a stepwise fashion (lowest contributor removed in each step until AUC stopped increasing). Bio15, Slope, land use, percentage of coniferous forest cover, percentage broadleaved forest cover, percentage of open water and density of roads were all removed during this process.

Maxent species distribution models
A Maxent model was performed in Maxent v 3.2.18. The last model consisted of 10 replications performed using Bio06 (minimum temperature of coldest month), Bio14 (mean precipitation of driest month), MaxTempChange, MinTemp-Change and Elevation. Jackknife tests and response curves were used to assess variable importance and responses to environmental factors.
The model was trained using records from all countries of E quercinus' known former range, except Russia (due to the lack of detailed data on the species' historical distribution in this region)-this also helped prevent overfitting of the model (Anderson and Raza 2010). The model was performed using with cross validation (where the presence points are divided into folds for training and test data for each replicate) and 10 replicate runs. 10,000 background points were used, due to earlier findings that such large background samples improve the accuracy of Maxent (Barbet-Massin et al. 2012).
All 10 models were assessed using the area under the receiver operating characteristic curve (AUC), with AUC values of 0.5 indicating the model is equivalent to a random model and values of 1 indicating a model with excellent fit for the data. The minimum standard set for this model was 0.7 AUC as suggested by Hijmans (2012).
Final variable importance was also assessed in the Maxent output. The percentage contribution was calculated by adding (or subtracting, if negative) the change in normalised accuracy gain to the contribution of the corresponding variable for each iteration of the model. For the calculation of permutation importance, for each environmental variable in turn, the values for the variables on training presence and background data are permuted at random and subsequently re-evaluated on the permuted data. The subsequent change in training AUC is then normalised to percentages. Jacknife tests, which display the improvement/reduction in model performance when the variable is omitted or used exclusively, were also performed. To create a binary map of E. quercinus distribution, a thresholded map was produced (using the calculated threshold value which maximised the sum of test sensitivity and specificity). Cells with suitability values above the threshold value were classed as suitable, whilst cells with values below it were classed as unsuitable.

Land use values per point
The "Extract Values to Points" tool in Arcmap 10.6 was used to find the values for land use [categories taken from the Corine land cover 2018 (CLC) 100 m 2 resolution raster (Copernicus Project 2020)] and HNVF status (true/false taken from the HNVF 2006 raster) (Paracchini and Oppermann 2011) for every presence point. Land use types were then simplified to improve analysis. CLC values for continuous/discontinuous urban fabric, industrial/commercial units, ports, airports and construction sites were all pooled under a single category ("urban"). All non-forested agricultural land use types were pooled together as "agricultural", excluding fruit tree plantations and agroforestry. Broad leaved, coniferous and mixed forests were pooled as "forest". Grasslands, heathlands and sclerophyllous vegetation were pooled as "grassland/shrub". Green urban areas and sports/leisure areas were pooled into green urban areas. Bare ground and sparsely vegetated areas were merged as "Bare/sparsely vegetated". All other land use types were pooled under "other".

Results
The Maxent models produced a mean AUC of 0.781, indicating the models fit well to the data. The suitability values indicated in the map of raw suitability (Fig. 1) and thresholded suitability (Fig. 2) also concurred well with observations about the current species population trends. Specifically, the map for the present indicates high suitability values for France, Spain and parts of Italy, where the species is apparently relatively stable (Bertolino 2017). It should be noted that Portugal, where data were very scarce, was also predicted to have high relative suitability. Patches of low suitability were predicted for Germany, with suitability decreasing into eastern Germany where the species is known to be in decline (Bertolino 2017). Areas further to the north (including the Scandinavian countries) and east (including Belarus, Ukraine, Romania and Russia) all had low relative suitability.

Variable importance and response curves
Both the assessment of mean variable contribution across the 10 replicates of the final model (Table 1) and mean jacknife evaluation (Fig. 3) indicated that the most important variable Fig. 1 Mean relative habitat suitability for E. quercinus across Europe in the present according to Maxent modelling. Suitability was significantly lower across eastern regions and higher to the south and west Fig. 2 Thresholded map for E. quercinus, based on the suitability value that maximises sensitivity and specificity of the Maxent model. Cells above the threshold were those with suitability values higher than the threshold which maximised test sensitivity and specificity (0.1659). Areas above the threshold were largely confined to the west and south of Europe, with smaller areas in central Europe in model construction was Bio06 (minimum temperature in coldest month), with DEM (elevation) being the second most important. The two climate change indicators had notable contributions (MaxTempChange scored 21.5% and MinTempChange scored 19.3%), and outperformed Bio14 (precipitation of driest month)-the jackknife and variable importance scores indicated that these two novel variables in combination with conventional variables, improved the model output.
Variable importance was assessed using jackknife testswith models being tested on variations run only with the variable and question, and on variations with all other variables except the variable in question. Models run exclusively with Bio14 performed worst of all, while its omission (running models with only the other variables) only marginally worsened performance. For all metrics, omitting Bio06 reduced gain the most (indicating it had the greatest effect). Bio06 also possessed the highest gain when used in isolation for plots (B) and (C) in Fig. 3 below, indicating a high amount of useful information. MinTempChange had the highest gain in isolation based on training gain, although had less effect on test gain and AUC based measures. Higher values for Bio06 were associated with higher suitability, although the relationships between the other variables and suitability were more complex (generally indicating decreases in suitability beyond certain points, with the exception of Max-TempChange where this trend was reversed). Suitability did not change much in response to different levels of Bio14 (precipitation in driest month), although higher values did slightly positively influence suitability. Elevation (DEM) improved suitability up to a point, although decreased for very high-elevation areas. Results for the climate change indicator variables were mixed. Areas which had moderate increases to minimum temperature (MinTempChange) saw improved suitability, but suitability dropped sharply beyond a certain point. For climate change to maximum temperature (MaxTempChange), low to moderate levels of climate change caused a sharp fall in habitat suitability, although higher degrees of change improved suitability. Figure 3 graphically illustrates these results.

Land use at observed points
The proportion of points varied noticeably between land use category and presence of HNVF (Fig. 4). Agricultural (1260 points), forestry (923 points) and grassland/shrub (502 points) areas possessed the most records for E. quercinus. HNVF practices appeared to be important for E. quercinus presence in agroforest ecosystems and grassland/shrub ecosystems, as well as 33% of all points in agricultural areas (Table 2).

Discussion
Our Maxent model produced an adequate area under receiver operating characteristic curve (AUC) value, and map predictions concurred with expert accounts of E. quercinus current range and population trends (especially in Germany, where low suitability values match with reports of the species population and range decline). The few records located in Finland were an exception and resided in low suitability zones; these populations may be relicts. Our model indicated that the areas of highest relative suitability were located in western and southern Europe. Some pockets of suitable habitat were found well to the east of E. quercinus current known range, which could be surveyed for relict populations or studied for possible reintroduction.
Our results should be considered cautiously, since modelling principles require that the species is in equilibrium with its environment (Elith et al. 2010;Franklin et al. 2013), which appears to be the case for E. quercinus only across parts of its western range. This concern was mitigated by the timescale used (all presence points from before year 2000 were excluded) and model reliability was further improved by the large sample size. Nonetheless, the findings of the models must also be considered in concert with qualitative analyses and relevant literature.
An unexpected difficulty we found, when analysing the available data, was the different prevalence of records between nations. With exceptions, records were readily available for nations in western Europe (especially Spain, France, Italy, Switzerland and Germany). For example,  (Valente 2017). A similar lack of E. quercinus records was observed in Belgium (despite many records along the French side of the France-Belgium border) and Austria (despite extensive records on the opposite sides of the Italy-Austria and Switzerland-Austria border areas). These differences could have been caused by variations in the availability biodiversity information between countries. Northern nations such as Denmark, Norway, Sweden and Finland have active biological recording communities and effective infrastructure for biological data recording, yet we only found two recent (post 1999) records in Denmark and Finland and no records in Norway and Sweden, which indicated that the scarcity of E. quercinus records in these nations was probably linked with species rarity or absence. In Estonia, records of the species stopped after 1990, and the species is considered likely extinct in this and neighbouring countries (Juškaitis 2018). Previous studies suggested that the collapse of the species in eastern and northern Europe was genuine (Bertolino 2017;Meinig and Buchner 2012;Juškaitis 2018), not an artefact of lack of sampling. An unfortunate limitation of our study was the lack of granular spatial data on historical land use practices. Thus, while we were able to calculate relative climate change across the study area, we were unable to do this for land cover. Although some historical land use data are available Fig. 3 Jacknife of training (A), test (B) and AUC (C) gain for each variable. The dark blue shows how test models performed (as measured by AUC) with only that specific variable (see variable name on left), while the teal blue indicates how the model performed without that variable from 1992 (European Space Agency 2020), E. quercinus' range collapse predates this by a large margin. Furthermore, remote-sensing-based identification of land use/cover (LUC) did not necessarily reflect land use practices within LUC categories (pesticide use and the removal of hedgerows being particularly relevant). For example, one "pixel" of agricultural land may have had different provision of wildlife habitats or pesticides in different time periods, or different provision than even neighbouring agricultural pixels in the present. The present landscape variables we modelled (land use category, density of roads, percentages of broadleaved and coniferous forest cover and slope) were all eliminated in early refinement runs due to their low importance. This cannot be considered evidence against land use change as a driver of E. quercinus range loss. Indeed, the relatively high proportion of presence points within HNVF areas could indicate the utility of nature-friendly farming practices for this species. In fact, the low scores (and subsequent elimination) of landscape variables may have been caused by either a real-world lack of importance relative to climate variables, or by the aforementioned limitations in the available GIS data (data deficiency regarding management practices and history per pixel).
Although the use of our novel variables (MinTemp-Change and MaxTempChange) offered a potentially interesting insight into E. quercinus biology and could be applied to other species in subsequent studies, we only used them to generate prediction of current habitat suitability. Calculating "future" versions of these variables was conceivable, by calculating the difference between present climate and the climate predicted by general circulation models. We chose not to in this case, both because our variables were novel and because non-climatic factors may also have driven changes in E. quercinus' range. To our knowledge, we were the first to incorporate historical recent climate change into an SDM. We argue that this approach could be of potential use in identifying species threatened by climate change. These novel variables could provide an intriguing alternative to dedicated eco-evolutionary methods, such as those proposed by .
Previous studies on other small mammals identified elevation shifts that could be related to climate change (Kelt and Meserve 2014). The observed variable response curves from our model showed that in areas of high increases in minimum temperature, suitability decreased sharply, while the opposite trend was observed for high increases in maximum temperature. One possible explanation for these non-linear  curves was that warming climates may have forced/allowed E. quercinus to colonise higher elevations, although this possibility would need to be verified by field studies. It is plausible that local populations of E. quercinus are adapted to locally specific climatic conditions and patterns (for example, their hibernation cycles could have become linked to locally specific averages of winter temperature). The declining E. quercinus population in southern Spain could indicate a vulnerability to local climate change (Santoro et al. 2017). Climate change has been observed to trigger physiological responses (Parmesan 2006) and or environmental/ trophic changes in other species (Gilman et al. 2010). Climate change has been observed to be especially disruptive during winter months (Williams et al., 2015). This could be particularly troublesome for E. quercinus whose fecundity and life history are affected by the onset of winter hibernation (Mahlert et al. 2018). If further evidence of climate change effects on E. quercinus is discovered, then reintroduction efforts for this species could use general circulation models to predict future climate patterns and thus identify ideal reintroduction areas. The hypothesis regarding climate change could be tested with longer term monitoring of specific populations. Remote methods to analyse the effect of climate on distribution patterns could be of significant use. If climate change is indeed a driver of range loss, then there is merit to attempting to model the species' response to predicted shifts in climate.
As a presence-background method, Maxent cannot truly estimate the probability of E. quercinus persisting in an area, instead it estimated the relative (not absolute) probability of observation/occurrence (due to imperfect detection and sampling effort bias). Our use of a target-group sampling bias file helped to significantly reduce the effect of differing sampling effort across the landscape (and therefore helped to reduce effect of imperfect detection, where detection is influenced by broader environmental trends), although even this was not a perfect solution because it only indirectly inferred sampling effort (Guillera-Arroita et al. 2015). The unknown prevalence (proportion of suitable sites occupied by the species) of E. quercinus rendered a calculation of true probability of occurrence/presence impossible. This issue further impacted our "thresholded" map ( Fig. 2) and the areas above threshold should therefore not be considered a perfect representation of the E. quercinus distribution. Nonetheless, our extensive efforts to correct for bias should mean our models produced a strong map of relative suitability and reliably showed approximate relationships with environmental variables. If true presence-absence data become available for even a small region, presence-absence modelling methods could further improve our understanding of E. quercinus distribution in the wider landscape. Similarly, ensemble presence-background modelling approaches with other algorithms could prove useful in future; especially as methods for reducing sampling bias improve and more data become available in currently low-sampled areas.
The high percentage of records observed within HNVF areas could be considered a promising indicator of the value of nature-friendly farming practices in small mammal conservation. Specifically, the high share of HNVF in the results for agroforestry (88%) and grassland ecosystems (74%) may indicate a role for nature-friendly agricultural areas in future conservation initiatives for this species. However, it cannot be ruled out that HNVF areas may have had a higher detection probability for this species, indeed biodiversity data were one of the criteria behind HNVF designation (Campedelli et al. 2018). The superior coverage of HNVF in the western European countries compared to most areas in eastern Europe could have played a role in the species continued persistence in these countries. However, data on HNVF are spatially limited and did not cover the whole Maxent training extent (which meant it could not be included in Maxent modelling), although all our presence points were within the HNVF data coverage. In addition, HNVF identification methodology relied on indirect methods (indicators such as important bird/butterfly areas and protected areas) rather than data on actual management practices. Certain management practices may also be more threatening to E. quercinus than others, for example, the removal of rocks (reducing available shelter) (Mori et al. 2020) and pesticide use (decreasing food availability).
It is possible that land use change, habitat degradation and fragmentation, management practices and climate change interacted to cause the decline witnessed in the species. By extension, favourable management strategies and improved habitat connectivity could help to reduce pressure on the species even if climate change is indeed a major driver behind E. quercinus decline. These measures, especially improved habitat connectivity, could help prevent further range losses by allowing time for evolutionary rescue (i.e., giving populations of E. quercinus time to adapt to new conditions through natural selection) (Diniz-Filho and Bini 2019).
In conclusion, our use of historical climate change data could prove useful in studies of other declining species. This study also provided some evidence suggesting that agricultural intensification and land use could be major drivers behind the decline of this species, and for the value of HNVF in conserving E. quercinus. Climate change and agricultural intensification could both play a role, or act synergistically, to explain the current distribution and status of the species. This study should not be seen as conclusive evidence; however, and subsequent field/monitoring studies should be performed to identify which environmental and landscape variables are the most important drivers.
Acknowledgements Dr. Rosalind Kennerly (IUCN SMSG, Bath, United Kingdom) for her advice on the possible causes of Eliomys quercinus decline. Kip D'arcourt and Paul Wexler (Wiltshire Mammal Group) for their insights on dormouse biology. We thank Flavia Tirelli and Maria Joao Ramos Pereira (BiMaLab, UFRGS, Brazil). for their advice on species distribution modelling. We thank Thiago Silveira (Federal University of Santa Catarina, Brazil) for his advice on modelling. We thank Victoria Graves for her suggestions to the final paper. We also thank the editor Frank Zachos, associate editor Mauro Lucherini and two anonymous reviewers for their feedback and constructive comments. The project was financially supported by the European Commission through the program, Erasmus Mundus Master Course-International Master in Applied Ecology (EMMC-IMAE) (FPA 2023-0224/532524-1-FR-2012-1-ERA MUNDUS-EMMC)-Coordinator F-J Richard, Université de Poitiers.
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.