Valentin Journé
email: valentin.journe@inra.fr

Jean-Yves Barnagaud

Cyril Bernard

Pierre-André Crochet

Xavier Morin

Pierre- André Crochet

Correlative climatic niche models predict real and virtual species distributions equally well

Keywords: Climate modelling of virtual species 4 Figures, 1 Table, 1 Supporting Information climate change, model evaluation, null models, range shift, spatial autocorrelation, species distribution models, blocked cross-validation

INTRODUCTION

Climate is one of the main drivers of global diversity patterns along latitudinal and elevation gradients [START_REF] Humboldt | Essai sur la géographie des plantes[END_REF][START_REF] Merriam | Laws of temperature control of the geographic distribution of the terrestrial animals and plants[END_REF][START_REF] Francis | A Globally Consistent Richness-Climate Relationship for Angiosperms[END_REF][START_REF] Willig | Latitudinal gradients of biodiversity: patterns, scale, and synthesis[END_REF]) and is unquestionably a major determinant of species distributions regardless of taxa (Woodward 1987, Kearney and[START_REF] Kearney | Mechanistic niche modeling: Combining physiological and spatial data to predict species' ranges[END_REF], as shown by paleoecological data (e.g. [START_REF] Williams | Late-Quaternay vegetation dynamics in North America: Scaling from taxa to biomes[END_REF], experimental approaches (e.g. [START_REF] Pigott | Factors controlling the distribution of Tilia cordata at the Northern limit of its geographical range. II. History in North-West England[END_REF]Huntley 1980, Rehfeldt et al. 2002) or current observations of range shifts in response to climate change (e.g. [START_REF] Walther | Ecological responses to recent climate change[END_REF][START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF]. However, a complex array of factors interacts with climate to shape actual species distributions. Biotic interactions, dispersal, micro-habitat availability, disturbances (including anthropogenic habitat changes or local species extirpation) all contribute to displace boundaries of species ranges away from the limits of their climatic niche [START_REF] Sexton | Evolution and ecology of species range limits[END_REF]) and generate differences between the climatic requirements of a species (i.e. its fundamental climatic niche in the sense of [START_REF] Hutchinson | Concluding remarks[END_REF]) and its geographic realization [START_REF] Gaston | The Structure and Dynamics of Geographic Ranges[END_REF].

Deciphering the links between climate and species distributions are more crucial than ever: robust predictions of species range shifts in response to climate change is a priority for the conservation of biodiversity and its functions [START_REF] Coreau | The rise of research on futures in ecology: rebalancing scenarios and predictions[END_REF][START_REF] Mouquet | Predictive ecology in a changing world[END_REF], but climate-based predictive models of species distributions rely on crucial assumptions on the role of climate in setting species distribution boundaries [START_REF] Sexton | Evolution and ecology of species range limits[END_REF]). These assumptions are straightforward for species whose range limits are known to be affected by climate, but they require a cautious examination in most cases because non-climatic drivers may be important [START_REF] Guo | Plant abundance: The measurement and relationship with seed size[END_REF]. In fact, disentangling the relative influence of climate and other ecological factors on species distribution, and understanding how these drivers determine species range limits, remains particularly challenging [START_REF] Coreau | The rise of research on futures in ecology: rebalancing scenarios and predictions[END_REF].

Various approaches have sought to identify the determinants of range limits in plants and animals [START_REF] Gaston | The Structure and Dynamics of Geographic Ranges[END_REF][START_REF] Sexton | Evolution and ecology of species range limits[END_REF]). Among these, frameworks based on correlations between species distributions and climatic variables have a long history [START_REF] Merriam | Laws of temperature control of the geographic distribution of the terrestrial animals and plants[END_REF]. Building on this legacy, correlative species distribution models (cSDMs, see [START_REF] Jarnevich | Caveats for correlative species distribution modeling[END_REF] represent the most popular approach by far and have been applied to a wide range of taxa and geographic scales (e.g. [START_REF] Thomas | Extinction risk from climate change[END_REF][START_REF] Pereira | Scenarios for Global Biodiversity in the 21st Century[END_REF][START_REF] Thuiller | Consequences of climate change on the tree of life in Europe[END_REF]. Experiments permit a more direct assessment of the causal role of climatic factors on species distributions (e.g. [START_REF] Breeman | Relative importance of temperature and other factors in determining geographic boundaries of seaweeds: Experimental and phenological evidence[END_REF][START_REF] Loehle | Height growth rate tradeoffs determine northern and southern range limits for trees[END_REF], Morin et al. 2007a), but they are difficult to carry out at large spatial scales and for all organisms. Mechanistic modeling approaches relying on explicit species responses (eg. growth, survival, reproduction, phenology) to climate variables have also been developed in the last two decades, for plant (e.g. [START_REF] Chuine | Phenology is a major determinant of tree species range[END_REF] and animal species (e.g. [START_REF] Kearney | Mechanistic niche modeling: Combining physiological and spatial data to predict species' ranges[END_REF], emphasizing the role of acclimation to cold temperatures [START_REF] Leinonen | A simulation model for the annual frost hardiness and freeze damage of Scots Pine[END_REF], phenology [START_REF] Chuine | A unified model for budburst of trees[END_REF], or metabolic rates [START_REF] Kearney | Mapping the fundamental niche: Physiology, climate, and the distribution of a nocturnal lizard[END_REF]. However, these process-based models have been restricted to a small number of species (see e.g. Morin et al. 2007b[START_REF] Sinervo | Erosion of lizard diversity by climate change and altered thermal niches[END_REF][START_REF] Diamond | A physiological trait-based approach to predicting the responses of species to expiremental climate warming[END_REF] because the data needed to implement such tools are difficult to collect.

The vast majority of studies aiming at predicting species range shifts in response to climate change have therefore used cSDMs based on climatic variables (Guisan andThuiller 2005, Pacifici et al. 2015). These "climatic cSDMs" predict local presence probabilities of a species from correlations between distribution data (occurrences in most cases, sometimes abundances) and climatic variables [START_REF] Elith | Species Distribution Models: Ecological explanation and prediction across space and time[END_REF] and are relatively straightforward to calibrate.

Non-climatic predictors can also be used in cSDMs (e.g. land cover: [START_REF] Thuiller | Do we need land-cover data to model species distributions in Europe?[END_REF][START_REF] Ay | The economics of land use reveals a selection bias in tree species distribution models[END_REF], dispersal and biotic interactions: Heikkinen et al. 2007[START_REF] Cabral | Estimating demographic models for the range dynamics of plant species[END_REF][START_REF] Boulangeat | Accounting for dispersal and biotic interactions to disentangle the drivers of species distributions and their abundances[END_REF][START_REF] Palacio | Biotic interactions in species distribution models enhance model performance and shed light on natural history of rare birds: a case study using the Straight-billed Reedhaunter (Limnoctites rectirostris)[END_REF], but these examples are the exception rather than the rule.

However, the ability of climatic cSDMs to properly identify and quantify the role of climate on species distributions is challenged increasingly often. The forecasting abilities of cSDMs can be tested by comparing their predictions of range shifts between past and current distributions -assessed with historical occurrence databased on observed climate change with actual range shifts. Although there are still few examples [START_REF] Elith | Species Distribution Models: Ecological explanation and prediction across space and time[END_REF], this approach yielded mixed results so far: [START_REF] Sofaer | Misleading prioritizations from modeling range shifts under climate change[END_REF] reported that modelling was poor at predicting observed range change in North American birds, while [START_REF] Fordham | How complex should models be ? Comparing correlative and mechanistic range dynamics models[END_REF] concluded the opposite in British birds (but with considerable variation in predictive abilities among species, see e.g. Fig. 3 in [START_REF] Fordham | How complex should models be ? Comparing correlative and mechanistic range dynamics models[END_REF]. For butterflies, [START_REF] Kharouba | Historically calibrated predictions of butterfly species' range shift using global change as a pseudo-experiment[END_REF] showed that climate changeinduced range shifts of North American species could be well predicted by SDMs for most of the 160 species tested but also reported a large variation between species in the abilities of the models to correctly predict observed range shifts despite high spatial model accuracies during model training. Another widely used approach to evaluate predicting abilities of cSDMs is to compare predicted and observed distribution in invasive areas. Again, this approach produced contrasted conclusions. For instance, Peterson et al. (2003, plants); Petitpierre et al. (2012, plants) or Ramírez-Albores et al. (2016, plants) reported good abilities of cSDMs in predicting invasive ranges but Broennimann et al. (2007, plants) ;Beaumont et al. (2009, plants); Gallagher et al. (2010, plants); Early and Sax (2014, plants); Goncalves et al. (2014, plants) or Hill et al. (2017, insects) found opposite conclusions.

One of the main reasons for the success of climatic cSDMs undoubtedly lies in the strong accuracy of their results when compared to observed species distributions. Such high goodnessof-fit under current climatic conditions is often interpreted as evidence for the ability of cSDMs to represent a causal link between climatic variables and species range (e.g. Pearson andDawson, 2003, 2004, and see review by [START_REF] Elith | Species Distribution Models: Ecological explanation and prediction across space and time[END_REF] or as a way to identify key specific climatic variables driving species range limits (e.g. [START_REF] Bombi | Modeling Bedriaga's rock lizard distribution in Sardinia: An ensemble approach[END_REF][START_REF] Gogol-Prokurat | Predicting habitat suitability for rare plants at local spatial scales using a species distribution model[END_REF][START_REF] Javanbakht | Genetic diversity and Quaternary range dynamics in Iranian and transcaucasian tortoises[END_REF]. However, several studies have questioned the ability of cSDMs to identify accurately the role of climate in determining species range [START_REF] Loehle | Model-based assessment of climate change effects on forests: a critical review[END_REF][START_REF] Bahn | Can niche-based distribution models outperform spatial interpolation?[END_REF][START_REF] Beale | Opening the climate envelope reveals no macroscale associations with climate in European birds[END_REF][START_REF] Lozier | Predicting the distribution of Sasquatch in western North America: Anything goes with ecological niche modeling[END_REF][START_REF] Fourcade | Paintings predict the distribution of species, or the challenge of selecting environmental predictors and evaluation statistics[END_REF][START_REF] Rich | Are North American bird species' geographic ranges mainly determined by climate?[END_REF].

A key criticism is that cSDMs rely on correlations that do not necessarily have a causal basis, and that they cannot be externally validated since they are typically tested by resampling the dataset of available occurrences [START_REF] Dormann | Correlation and process in species distribution models: Bridging a dichotomy[END_REF].

Null models represent a more appropriate way to assess the capacity of species distribution models to detect the causal relationship between species and climate [START_REF] Gotelli | Null models in ecology[END_REF].

Only a handful of studies have used null models to assess the robustness of climatic cSDMs in the past, either by comparing the output of cSDMs obtained on actual species distributions and on virtual species distributions [START_REF] Bahn | Can niche-based distribution models outperform spatial interpolation?[END_REF][START_REF] Beale | Opening the climate envelope reveals no macroscale associations with climate in European birds[END_REF], Jiménez-Valverde et al. 2011) or by using randomized environmental predictors (Boucher-Lalonde and Currie 2016, [START_REF] Fourcade | Paintings predict the distribution of species, or the challenge of selecting environmental predictors and evaluation statistics[END_REF]). These null model approaches have yielded mixed results. For instance, [START_REF] Beale | Opening the climate envelope reveals no macroscale associations with climate in European birds[END_REF] concluded that species-climate associations found by cSDMs were no better than chance in a majority of European bird species, while [START_REF] Jiménez-Valverde | Dominant climate influences on North American bird distributions[END_REF] found the opposite conclusion for North American birds. However, most of these studies questioning the actual robustness of cSDMs have focused on bird species, and no study that we know of has used null models to assess the goodness-of-fit of cSDMs across several taxonomic groups, i.e. to test whether the robustness of cSDMs may depend on the group of taxa considered.

In this study, we used real species from five taxonomic groups and virtual species distributions of the same range size to test whether or not the accuracy of cSDMs calibrated with climatic variables reflect the role of climate in shaping species distribution. By doing so, we assessed the extent to which cSDMs can actually capture the processes linking climate and species distributions. A key requirement of our approach is the design of appropriate null models. We used random diffusion models, unconstrained by climatic gradients, to build virtual species distributions that have no causal link to climate, then compared the performance of climatic cSDMs when predicting real and virtual species distributions. We selected real species in five groups differing in their ecology and dispersal abilities, which excludes that our findings are explained by biological traits specific to a given group of organisms. Our aim was to answer the following question: can the distribution of virtual species be predicted with a similar accuracy as real species? Doing so, we do not question the role of climate as a major driver of species distributions (see [START_REF] Merriam | Laws of temperature control of the geographic distribution of the terrestrial animals and plants[END_REF][START_REF] Walther | An ecological "footprint" of climate change[END_REF][START_REF] Chen | Rapid Range Shifts of Species[END_REF], or that cSDMs can be useful for assessing climate change impacts whenever climate is the main factor limiting a species range. Instead, what we precisely question is the use of cSDMs to assess the importance of climate on species distributions and consequently their use in predicting range shifts in response to climate change in species where the link between climate and range limits has not been clearly established.

We focus this study on the Western Palearctic region (Europe, North Africa and the Middle East) because species distributions there are known with greater accuracy than in most other areas and can be associated with reliable climatic data. We selected two plant functional groups (herbaceous plants and trees), and three vertebrate groups, two ectotherms (reptiles and amphibians) and an endotherm (birds). We randomly picked 25 species (hereafter 'real species') of widely variable range size (i.e. prevalence) in each of these five groups. For each real species, we generated 10 virtual species distributions with spatial cohesion. We hypothesized that the predictive power of cSDMs varies across each taxonomic group due to their contrasted ecologies.

For example, one may expect that plants, amphibians and reptiles could be more influenced by the variation in climate variables than warm-blooded species such as birds [START_REF] Gaston | The Structure and Dynamics of Geographic Ranges[END_REF]. We also expect that the predictive accuracy of cSDMs decreases with increasing species range size, because larger distributions tend to sample broader arrays of climatic conditions than restrictedrange species. Furthermore, we used a newly-proposed rigorous evaluation of models' outcomes [START_REF] Roberts | Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure[END_REF] to strengthen the scope of our findings. We therefore i) assessed the sensitivity of our analysis to range size and ii) compared accuracy of cSDM results on real species distributions between the five species groups, then iii) compared the accuracy of climatic cSDM suitability maps obtained for these virtual species distributions to those obtained for the real species, across the five groups of taxa.

METHODS

Study area and climatic variables

The study area corresponds to the Western Palearctic ecozone, including the European region of Eurasia together with North Africa (north of the Sahara) and the Middle East (Fig. 1). It is located roughly between the latitudes of 30° to 72° N and the longitudes of 12° W to 68° E and covers an area of 33,152,400 km² (i.e. 331,524 S1) and were retained as predictors in all subsequent modeling analyses.

This was sufficient for our purpose as we did not aim to identify the most relevant climatic variables influencing the range of a given species.

Real species distribution data

We first built a species list for birds, amphibians, reptiles and herbaceous plants occurring in Europe, retaining only species native to Europe that have most of their breeding range within Europe. We further excluded seabirds and other strictly coastal or insular bird species. Once the lists were completed, we randomly selected 25 species for each of these taxa (100 species in total).

For these four taxonomic groups, species distributions were downloaded from the IUCN Red List website (http://www.iucnredlist.org). If a distribution was not available for a given species, we randomly picked another species of the same group in our list. For trees, we downloaded the distributions of all 26 available species (at the time of the study in 2016) from EUFORGEN (http://www.euforgen.org, San-Miguel-Ayanz et al. 2016). We removed areas of non-native presence (if any) and non-breeding range for birds, as well as parts of the range located outside the study area. We then rasterized species distribution maps into presence-absence data at a 5minute resolution to match climatic data. Because some of the selected species had small ranges (i.e. less than 100 pixels), we anticipated difficulties to fit accurately a distribution model (as illustrated in other studies, e.g. [START_REF] Stockwell | Effects of sample size on accuracy of species distribution models[END_REF]. We therefore added three supplementary species of amphibians and herbaceous plants to circumvent this issue. In total, we gathered data for 132 species distributions of five taxonomic groups. The distributions of 88 species were entirely included within our study zone. The distributions retained for the final analyses vary from local endemics (minimum of six pixels) to continent-wide presence (maximum of 244,492 pixels on 331,524 land pixels, see Metadata S1).

Null models

For each real species, we generated 10 virtual species distributions with the same range size (prevalence) as the real species (Fig. S1, Fig. 1). We chose to generate spatially cohesive distributions for virtual species, i.e. the virtual distributions were formed of a single block rather than multiple isolated pixels with occurrence data (see discussion). We developed an algorithm under the R environment (R Core Team 2016) to generate these virtual species distributions using a diffusion model approach (Metadata S2). We started by randomly selecting one pixel in the study area. The virtual distribution was then extended gradually by colonizing adjacent pixels in a direction chosen randomly at each step, until the final range size was reached (i.e. the same range size as the corresponding real species). Whenever a pixel was selected within a large water body, it was exchanged with the closest unoccupied pixel in land. This procedure guaranteed that sea basins were not a barrier to "colonization" and avoided the mid-domain effect [START_REF] Colwell | The mid-domain effect: Geometric constraints on the geography of species richness[END_REF], where most distributions would be located in the center of the study area.

Distribution modeling

All cSDMs were fitted with the R package BIOMOD2 [START_REF] Thuiller | BIOMOD -A platform for ensemble forecasting of species distributions[END_REF]). We used the same modeling workflow for real and virtual distributions (Fig. S1). We used two kinds of algorithms to fit the cSDMs. First we built generalized additive models (GAMs; [START_REF] Wood | Generalized additive models: an introduction with R[END_REF] to describe the relationships between presence-absence data and the six climatic predictors. GAMs are commonly used in species distribution modeling because they handle complex responses and are robust to data quality as well as to the number of absence data [START_REF] Barbet-Massin | Selecting pseudo-absences for species distribution models: How, where and how many?[END_REF]).

As our main aim was to compare real with virtual species, we did not need to split the dataset of observed presences as commonly done to test for the validity of the fitted model. Thus we first fitted GAMs by including all the pixels in our study area (331,524 pixels, from six to 244,492 presences depending on the species), considering all pixels outside species' distributions as absences.

However, considering all pixels may inflate goodness-of-fit values and thus level differences between real and virtual species. In most cSDMs-based studies, model validation is assessed on a random subsample of the dataset of observed presences or presences-absences different from the subsample of data used for model calibration. Calibration and validation based on random samples of presence-absence data can underestimate the errors when dependence structures are present in the data [START_REF] Roberts | Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure[END_REF]), as it is the case with climatic data. In fact, defining separate blocks

for training and testing model represents a new opportunity to limit effect of spatial structure for climate and species distribution [START_REF] Roberts | Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure[END_REF]. Block cross-validation can address these issues. As recommended by [START_REF] Roberts | Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure[END_REF], we defined two non-randomly chosen blocks of data, one for cSDM calibration and another one for cSDM validation. To do so, we split the study area in two according to the mean latitude of species distribution. We then sampled 5000 absences in each block to insure a minimal value of 10,000 absences in total for the study site [START_REF] Barbet-Massin | Selecting pseudo-absences for species distribution models: How, where and how many?[END_REF]. When it was impossible to split the study area according to latitude, we used longitude. We present results for both methods, i.e. with and without cross validation, for GAM algorithm. Our general methodology is summarized in Figure S1.

As GAMs are prone to overfitting and may produce response functions that are difficult to interpret biologically [START_REF] Wood | Generalized additive models: an introduction with R[END_REF], we also used MAXENT (v3.3.3, Phillips et al. 2006, using only presence data) to assess the robustness of our findings to changes in algorithm. The results obtained with MAXENT were qualitatively similar to those obtained with GAMs and are not detailed here (see Appendix A1 in Supporting Information). Note that overfitting can also affects machine learning methods or MAXENT [START_REF] Elith | Novel methods improve prediction of species' distributions from occurrence data[END_REF][START_REF] Elith | Do they? How do they? Why do they differ? on finding reasons for differing performances of species distribution models[END_REF].

Models evaluation and statistical analysis

We used the two most common evaluation methods [START_REF] Fourcade | Paintings predict the distribution of species, or the challenge of selecting environmental predictors and evaluation statistics[END_REF]) to calibrate and validate GAMs, the Area Under the Receiver Operating Curve (AUC) criterion and the TSS. AUC is independent from the threshold value used to convert probabilities of presence per pixel into presenceabsence data [START_REF] Elith | Species Distribution Models: Ecological explanation and prediction across space and time[END_REF] and is one of the most popular criteria to evaluate cSDMs. An AUC value of 1 means a perfect fit, while a value of 0.5 corresponds to random prediction of a species presence in a given cell.

To better check whether our conclusions are sensitive to the evaluation metrics used, we also used a threshold-dependent metric with the GAM algorithm, the True Skill Statistics (TSS, [START_REF] Allouche | Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic (TSS)[END_REF]). The TSS is an alternative metric used in cSDMs studies that is not influenced by the distribution size. It is based on presence-absence data and varies between -1 (poor fit) and 1 (perfect fit). To convert the probabilities of presence into presence-absence data for TSS, we used a threshold value maximizing sensibility and specificity as recommended by [START_REF] Liu | Selecting thresholds for the prediction of species occurrence with presence-only data[END_REF].

Two real species and 49 virtual species were not predicted due to the low number of presences used to calibrate the model, but the results yielded the same conclusions than the AUC (see Figure 3).

To compare the accuracy of cSDMs results for virtual and real species distributions, we assigned a rank to each real species based on its AUC (or TSS) value compared with the AUC (or TSS) values of its 10 corresponding virtual species: if the real species was found to have the highest AUC (or TSS) it was assigned the first rank, i.e. rank 1, while if it was found to have the lowest AUC (or TSS) it was assigned the last rank, i.e. rank 11. We then compared the rank distribution of the real species to the expected distribution under the null hypothesis that AUC (or TSS) values did not differ between real and virtual species. We also used a t-test to test whether AUC (or TSS) values of cSDMs differed between real and virtual species and an F-test of equality of variances I to test if real and virtual species had AUC (or TSS) of equal variance (Stats package in R).

Last, we ran a General Linear Model in Statistica version 10 (StatSoft 2011) to examine the effects of group and prevalence on AUC (or TSS) for real species only. Prevalence was log-transformed prior to analysis to reduce non-normality and group was defined as random. We started with the complete model AUC (or TSS) = log(prevalence) + group + group * prevalence.

RESULTS

Comparison between real and virtual species

The distributions of the real species were generally accurately predicted without cross-validation (hereafter "CV") with AUC values ranging from 0.899 to 1.00 (µAUC = 0.972, σAUC = 0.02), while AUC values with blocked CV ranged from 0.308 to 0.99 (µAUC-CV = 0.675, σAUC-CV = 0.179, Figure S3). The TSS values ranged from 0 to 0.99 without CV (µTSS = 0.844, σTSS = 0.129) and from 0 to 0.98 with blocked CV (µTSS-CV = 0.323, σTSS-CV = 0.278, Figure S3). The AUC and TSS values of cSDMs for virtual species distribution were also higher without CV (µAUC = 0.975, σAUC = 0.05 ; µTSS = 0.848, σTSS = 0.20), than with blocked CV (µAUC-CV = 0. Most real species distributions had a lower rank than virtual species distributions, and only eight of the 132 real species were ranked first in comparison to their virtual species. The majority of virtual species distributions were thus predicted with a stronger accuracy than real species distributions (Fig. 2, Fig. 3-a & Metadata S1). The proportion of real species ranked first is very similar to the expected proportion if species rank was randomly assigned, regardless of their "status" (i.e. real or virtual): under random expectations, one out of 11 real species should show a better fit than its 10 associated virtual species, i.e. a total of 12 (132/11) real species ranked first by chance. The results were similar regardless the metric and the methodology used (Fig. 3-b, c, d & Metadata S1) or the algorithm (Figure S2 & Metadata S1). We can thus conclude that real and virtual species' distribution were predicted with the same accuracy.

Effect of prevalence and taxon on the fit of the cSDMs for the real species

Prevalence strongly and negatively affected the cSDMs accuracy (AUC and TSS, without CV and with blocked CV) for real species (Table 1). The group did not have a significant effect on AUC and on TSS with blocked CV, but it significantly affected TSS without CV (Table 1). In this last case, the interaction between prevalence and group was also significant. In all cases, AUC or TSS seemed to decrease more for birds than for the other groups when prevalence increased (Fig. 4).

DISCUSSION

Our main result is that SDMs based on climatic variables (cSDMs) predict equally well the distribution of real species and of virtual species whose distribution is independent of climate, regardless of the method used to fit cSDMs. Using a blocked cross-validation showed that virtual distributions may be even better predicted with cSDMs than real species distributions. Our findings therefore demonstrate that a strong predictive power of cSDMs should not necessarily be interpreted as an accurate depiction of the link between climate and species distribution, and thus calls for a greater caution when using cSDMs to make predictions about future changes in species distributions in response to climate change.

Prevalence and taxonomic groups did not affect our conclusions

We expected the increase size of the distribution (or prevalence, see [START_REF] Stockwell | Effects of sample size on accuracy of species distribution models[END_REF]Peterson 2002, Proosdij et al. 2016) increase the accuracy of cSDMs, as previously found in invertebrates (Aguirre- [START_REF] Aguirre-Gutiérrez | Fit-for-Purpose: Species Distribution Model performance depends on evaluation criteria -Dutch hoverflies as a case study[END_REF]) and tetrapods (mammals, reptiles, amphibians and birds; Morán-Ordóñez et al. 2017). In agreement with these studies, we found that predictive accuracy of cSDMs decreased with increasing range size in real species. Contrary to our expectations, we did not find a significant variation of predictive performance of cSDMs between taxonomic groups, despite a tendency for weaker accuracies for birds (different slope of the AUCprevalence relationships, Fig. 4). We might speculate that this tendency for a lower accuracy of climatic cSDMs in birds would become significant with a larger sample size per taxonomic group (approx.

25 real species per group in our study), as expected if birds, being endotherms, are less constrained by climate than ectotherms and plants [START_REF] Gaston | The Structure and Dynamics of Geographic Ranges[END_REF]. The main objective of including several taxonomic groups in our study was not to test for differences in accuracy among groups, however, but to ensure that our main conclusion was not dependent on the taxonomic group used. Therefore, virtual species distributions can be predicted as accurately as real species distributions, and this pattern does not depend on the size of the distribution group and also appears consistent regardless the ecological characteristics of the real species investigated.

Spatial auto-correlation of climatic variables is the most likely culprit to explain our finding

How to explain that climatic variables can accurately predict a distribution range even when it is not linked to climate? [START_REF] Bahn | Can niche-based distribution models outperform spatial interpolation?[END_REF]

Deciphering the role of climate on species distributions

Most studies that used null models to assess the ability of cSDMS to identify causal links between climate and species distributions concluded that climatic cSDMs did not perform generally better than alternative null models (Bahn andMcGill 2007, Beale et al. 2008 Testing the role of climate on species distributions is more crucial than ever. If fitting cSDMs cannot inform us on the causal link between climate and species occurrence, which possibilities can be envisaged? Possible alternative approaches include mechanistic models linking ecophysiological and/or demographic components to climatic variables on the basis of empirical and/or experimental data [START_REF] Chuine | Phenology is a major determinant of tree species range[END_REF][START_REF] Sinervo | Erosion of lizard diversity by climate change and altered thermal niches[END_REF][START_REF] Gutiérrez | Using a dynamic forest model to predict tree species distributions[END_REF].

Although some of the processes embedded in these mechanistic approaches may be sometimes modelled partly through correlations (for instance if a process relies on bioclimatic thresholds inferred from presence-absence observed in distribution maps), the bias is certainly much weaker than for cSDMs that fundamentally depend on the species range size to be calibrated. Therefore, because they mostly rely on response functions that are a priori defined in experimental settings, they are not biased by spatial autocorrelation, and the effect of climate is directly tested. Another promising approach is the experimental manipulation of range limits (reviewed in [START_REF] Hargreaves | Are Species' Range Limits Simply Niche Limits Writ Large? A Review of Transplant Experiments beyond the Range[END_REF]. Last, the current level of climate change offers pseudo-experimental settings to directly test the impact of climate on species distribution [START_REF] Kerr | The macroecological contribution to global change solutions[END_REF]. Whenever actual changes in range limits as a response to climate change match the predictions of climatic cSDMs, the causal role of climate in determining distribution becomes well supported. As mentioned in the introduction, such approaches have yielded mixed responses within a general pattern of poleward and upward range shifts and it is likely that the direct contribution of climate to range limits vary between taxonomic groups and/or biogeographic regions (see [START_REF] Lenoir | Climate-related range shifts -a global multidimensional synthesis and new research directions[END_REF].

Unfortunately, none of these alternative approaches is as easy to use or as general as correlative SDMs.

Comparisons with previous studies and the need for repeatability in ecology

Our work is not the first one to use virtual species distributions to assess the power of cSDMs to capture actual links between climate and species distributions (see references in the introduction).

However, we wish to point out here the novel aspects of our study.

Firstly, the most original aspect of this work is clearly the comparison of several taxonomic groups of real species. Most previous papers using similar approaches focused on only one group of species (in most cases avian species), and it could be argued that their results were driven by a possible lack of response of birds to climate (in comparison with other kinds of organisms with physiological requirements more directly related to climate, like plants or ectotherms). In fact, one of the previous similar studies [START_REF] Beale | Opening the climate envelope reveals no macroscale associations with climate in European birds[END_REF] concluded that « the distributions of most birds in our study are not strongly associated with the climate variables currently available » rather than questioning the bases of correlative climatic niche modelling. We are not aware of any study comparing the performance of cSDMs on virtual and real species for such a wide range of taxa (ie., the four main groups of terrestrial vertebrates and two functional groups of plants), thus comparing species strongly differing in dispersal abilities (from birds and trees to salamanders) and ecophysiology (endotherms vs. ectotherms and plants). We detected some differences between groups (consistent with the ecology of the groups), yet our conclusions draw a strong general trend. Therefore, with studies such as ours, it will be difficult to argue that the lack of power of cSDMs is restricted to the group(s) investigated.

Secondly, validation practices rapidly evolve in niche modelling and low performance of cSDMs compared to null models has also been linked to failure of validation tools. The study by [START_REF] Fourcade | Paintings predict the distribution of species, or the challenge of selecting environmental predictors and evaluation statistics[END_REF] concluded on "the inability of current evaluation metrics to assess the biological significance of distribution models". Here, we used the most stringent and recent validation tools -in addition to the commonly used AUC -and we showed that our conclusions did not change with the validation tools.

Last, we also wish to point out that in spite of the previous studies mentioned above, common practice has not changed much in the field and cSDMs are still widely used for forecasting without much caution: cSDMs are more widely used than ever, and the limits of this modelling technique and philosophy are almost never acknowledged. Like in many other fields, it will take at least a series of independent studies by several teams (see [START_REF] Palmer | Quasi-replication and the contract of error: lessons from sex ratios, heritabilities and fluctuating asymmetry[END_REF][START_REF] Kelly | Replicating empirical research in behavioral ecology: How and why it should be done but rarely ever is[END_REF][START_REF] Nagakawa | Replicating research in ecology and evolution: feasibility, incentives, and the cost-benefit conundrum[END_REF] Parker 2015 on the need to repeat studies in ecology and evolution) before the cautionary message is heard. We thus think that it is important to publish independent studies reinforcing previous conclusions with different approaches and independent data sets. In that sense, we believe that our study based on several groups of species, relying on spatially-explicit virtual distributions (instead of randomizations) across Europe and using the most recent methods in cSDMs validation represents a significant step in this direction.

Methodological issues and perspectives of our study

The methods used in this study differ from the classical way of fitting cSDMs in several aspects.

Firstly, we used polygons of distribution range as presence data, resulting in the spatial segregation of presence and pseudo-absence pixels. Many studies using cSDMs used actual occurrence data instead, resulting in a mixture of presence and pseudo-absence data inside the distribution range of species. Such an approach assumes that i) pseudo-absence data inside the distribution range of species are more likely to correspond to real absences than to pixels where the species is actually present but has not been detected and ii) that the environmental factors used to model species occurrence differ between presence and pseudo-absence pixels. The first assumption might be realized in some species with high detectability in well-prospected areas but is certainly not verified in most empirical studies that we know of. The second assumption seems unrealistic in climatic cSDMs based on broad-scale climatic variables (such as the WORLDCLIM database), as most of the variables used for modelling certainly do not vary at the spatial scale corresponding to the spatial resolution of the presence and pseudo-absence data. In other words, it is highly unrealistic to use climatic variables varying over regional or continental scales to explain whether a species is present or absent in particular pixels of a local region. We are therefore confident that using polygons of distribution range as presence data does not undermine our work. We nevertheless checked the sensitivity of our results to this methodological choice by modeling real and virtual species distribution after removing 20%, 40%, 60% and 80% of the presence data and selecting pseudo-absence randomly inside and outside the species range, as done in many empirical studies (e.g. [START_REF] Araújo | Validation of species-climate impact models under climate change[END_REF][START_REF] Barbet-Massin | Selecting pseudo-absences for species distribution models: How, where and how many?[END_REF][START_REF] Elith | Do they? How do they? Why do they differ? on finding reasons for differing performances of species distribution models[END_REF]. Using various proportion of presence data did not alter our results and our conclusion that real species distributions are not predicted with stronger accuracy than virtual species distributions therefore remains valid (see Appendix A1 in Supporting Information).

Secondly, the size of the study area is known to affect the accuracy of modelling results of cSDMs [START_REF] Thuiller | Niche properties and geographical extent as predictors of species sensitivity to climate change[END_REF]. If the study area is too large relative to the species distribution, there is a risk of overfitting (i.e. increased AUC values) leading to less realistic estimates of niche characteristics (Anderson andRaza 2010, Barve et al. 2011). This is simply because increasing the range of climate conditions spanned in a study area increases the ability of cSDMs to discriminate the range of conditions suitable for a given species, but at the cost of decreasing the accuracy. In this work, we did not adjust the size of the study area to be proportional to each distribution (Fig. 1, Metadata S1) because defining a specific study area for every real and virtual species distribution would distort the modeling procedure and would not allow a fair comparison for all species. However, the key point was that real and associated virtual species always had the same range size. Furthermore, many real species selected had a large prevalence, so the effect of large study area relative to species range cannot be responsible for our main conclusion.

Thirdly, we chose to assess the performance of the model through a blocked crossvalidation [START_REF] Roberts | Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure[END_REF][START_REF] Fourcade | Paintings predict the distribution of species, or the challenge of selecting environmental predictors and evaluation statistics[END_REF]. This method provides a better spatial independence between calibration and validation datasets than usually done in cSDMs-based studies with random cross-validation. Doing so also strengthened the robustness of our results.

Lastly, the generated virtual distributions were on average less fragmented than the real species distributions (Fig. 1). This pattern arose because our null models relied on continuous distributions obtained through a diffusion model, while real species may have disjointed distributions, especially those located in mountains and/or at high latitude. It has been shown that the aggregation level may influence the goodness-of-fit of cSDMs, with species distribution with a high level of aggregation usually showing higher AUC values [START_REF] Beale | European bird distributions still show few climate associations[END_REF]. Although this approach is likely to have a weak impact on our results, a relevant perspective to this work would be to create a diffusion model allowing to design fragmented distributions, in order to compare real and virtual species' distributions with a similar level of spatial fragmentation and/or convolutions of distribution limits. However, we do not expect this to affect our main conclusion, because our set of real species also includes distributions that have simple shapes and they were not better predicted than virtual species distributions (for examples see Figure S4), and because, alternatively, some of the virtual distributions with a convoluted outline were adequately predicted by cSDMs. We are thus confident that the main finding of the present study is robust to this limitation about the level of fragmentation of the virtual distributions.

Combining various approaches to better predict the impact of climate change [START_REF] Thuiller | Do we need land-cover data to model species distributions in Europe?[END_REF][START_REF] Cabral | Estimating demographic models for the range dynamics of plant species[END_REF][START_REF] Boulangeat | Accounting for dispersal and biotic interactions to disentangle the drivers of species distributions and their abundances[END_REF]). Such efforts have improved the fit of the models to current distributions and shed light on the role of non-climatic drivers of species distributions, so they certainly represent a promising venue for forecasting the impact of climate change. The climatic part of these hybrid SDMs remain the same than in classic cSDMs in most cases, however, and their forecasting abilities should be rigorously evaluated, for example using historical distribution data or out-of-range introductions, as has been done with classical cSDMs.

Given these limitations, how can we obtain robust predictions of the impact of climate change on species distributions? As written above, although the use of process-based or experimental approaches is currently expending, they are not straightforward to develop. A first attempt to solve this issue may be to limit forecasting with cSDMs to species for which range shifts have been shown to be predicted with strong accuracy through calibration on historical distributions and validation with current distributions [START_REF] Kerr | The macroecological contribution to global change solutions[END_REF] or to invasive species whose invasive areas have been well predicted -although both approaches brought contrasted results as reviewed in the introduction. Other possibilities may be to focus on species where the link between climate and the probability of presence is grounded on solid biological bases, such as experiments in controlled conditions (Rehfeldt et al. 2002, Chuine and[START_REF] Chuine | Process-based models of phenology for plants and animals[END_REF]

TABLES

Table 1. Linear Mixed Model for the effects of group (random) and prevalence on evaluation metrics for GAM models with and without blocked cross-validation. The interaction term was removed from the final model when it was non-significant (p-values between 0.25 and 0.82 for the three models with non-significant interaction term). See text for detail.

CAPTIONS AND FIGURES

  or using measures of the effects of climate on life-history traits[START_REF] Pigott | Factors controlling the distribution of Tilia cordata at the Northern limit of its geographical range. II. History in North-West England[END_REF] Huntley 1980, Sinervo et al. 2010), and to select those climatic variables that are biologically meaningful prior to the modelling steps of cSDMs.Yee-Law et al. (2016) nicely illustrates how climatic cSDMs can provide accurate information on niche limits and suitability in a given location when experimental data demonstrate a causal link between climate and fitness; this was the case for 31 of the 40 plant species investigated, so for nearly a quarter (9 out of 40) of the species investigated a link between suitability as predicted by cSDMs and actual local persistence as assessed by transplant experiments was not supported. An even more mixed pattern emerges from studies using humanmediated species translocations as quasi-experiments (see Introduction section), so a direct link between climate and current range limits cannot always be taken for granted.CONCLUSIONDisentangling the role of climate and other factors in shaping species ranges overreaches fundamental ecology and has far-reaching consequences when trying to predict biodiversity responses to climate change. Fitting cSDMs using climatic variables and using them to predict range change in response to climate change rests on the untested hypothesis that climate determines where a species is found. Our study highlights that fitting correlative SDMs based on climatic variables to current species distributions is not enough to assess the effect of climate on species distribution. This especially calls for caution when forecasting the impact of climate change on species range. While cSDMs remain an important tool in the emergence of predictive ecology[START_REF] Mouquet | Predictive ecology in a changing world[END_REF], a better understanding of the drivers of species distribution is crucial for conservation and climate change impact issues. It will allow the development of more realistic forecasting models incorporating dispersal, biotic interactions, actual relationships between climate and fitness, and non-climatic environmental factors.

Figure 1 :

 1 Figure 1: Distribution of Chloris chloris (a, bird), Abies alba (b, tree) and Eryngium alpinum (c,

Figure 2 :

 2 Figure 2: Comparison of the AUC values without cross-validation (a) and with blocked cross-

Figure 3 :

 3 Figure 3: Ranks of metrics values (AUC & TSS) of all real species distributions modeled with

Figure 4 :

 4 Figure 4: Effect of taxonomic group and prevalence on AUC values.

  741, σAUC-CV = 0.16 ; µTSS-

CV = 0.432, σTSS-CV = 0.28). The goodness-of-fit values of models fitted without CV were similar between virtual and real species

(GAM AUC, df = 283.62, p = 0.206; GAM TSS, df = 200.16, p = 0.691)

. However mean values were different between real and virtual species distribution fitted with blocked CV (GAM AUC-CV, t-test, t = -4.36, df = 1450, p < 0.0001; GAM TSS-CV, t-test, t = -4.20, df = 1450, p < 0.0001), with virtual species showing higher AUC and TSS values on average.

The AUC and TSS values of cSDMs for virtual species distribution were more variable than for real species (AUC: F-test of equality of variances; F131, 1319 = 4.88, p < 0.00001; TSS: F-test of equality of variances; F131, 1319 = 2.39, p < 0.00001). However, these differences disappeared with blocked CV (AUC-CV: F-test of equality of variances; F131, 1319 = 0.832, p = 0.14; TSS-CV, Ftest of equality of variances; F131, 1319 = 1.04, p = 0.75).

  , Boucher-Lalonde and Currie 2016,Fourcade et al. 2018, Rich and[START_REF] Rich | Are North American bird species' geographic ranges mainly determined by climate?[END_REF] but see Jiménez-Valverde et al. 2011). Among these studies, some have argued that, because cSDMs based on climatic variables can predict null distributions with no relationship to climate as accurately as real species' distributions, or because environmental variables added no predictive power beyond what spatial interpolation could provide, climate might not be the main factor driving species' distribution (e.g.Beale et al. 2008, Rich and[START_REF] Rich | Are North American bird species' geographic ranges mainly determined by climate?[END_REF]. We do not share this view.In our opinion, the ability of climatic cSDMs to model virtual species' distributions with the same accuracy as real species' distribution tells nothing about the role of climate in shaping species' distributions. These results simply demonstrate that fitting climatic cSDMs is not necessarily an appropriate tool to test for the role of climate in shaping species distributions. Because climatic cSDMs can and will accurately predict distributions that have no link to climate(Bahn and McGill

2007,

Beale et al. 2008, our results)

, we need other approaches to investigate which species distributions are limited by climate or by other biotic or abiotic factors.

  Forecasting distributions (e.g. projecting climate change impact on species distribution) is conditioned by two assumptions: that current distributions are mainly limited by climatic variables and that climatic niches remain stable through time (this is called the "climatically constrained hypothesis" by[START_REF] Rich | Are North American bird species' geographic ranges mainly determined by climate?[END_REF]. It is probably sensible to assume that evolutionary changes in climatic niche will be rare enough (i.e. be expressed in few species) or small enough (i.e. be of low magnitude) at contemporary time scales to have a limited effect on most forecasting attempts, even if contemporary adaptation to new climate has been documented before (e.g.[START_REF] Colautti | Rapid Adaptation to Climate Facilitates Range Expansion of an Invasive Plant[END_REF][START_REF] Geerts | Rapid evolution of thermal tolerance in the water flea Daphnia[END_REF][START_REF] While | Adaptive responses to cool climate promotes persistence of a non-native lizard[END_REF]. However, assuming that all current species distributions are mainly limited by climatic variables is highly questionable, as explained in the introduction and earlier in the discussion. Our results further emphasize that forecasting distributions with cSDMs should be done with caution: correlative species distribution models may be useful to predict changes in species distribution due to climate change whenever climate is the main factor determining range limits[START_REF] Lee-Yaw | A synthesis of transplant experiments and ecological niche models suggests that range limits are often niche limits[END_REF]. For instanceKharouba et al.

	(2009) constructed cSDMs under historical conditions for butterfly species and predicted range
	shifts reasonably well when climate changed for many (but not all) species. This is consistent with
	the fact that the physiology of many insect species (phenology, survival rates) is strongly
	dependent on climate (Chuine and Régnière 2017). Since current distributions are not necessarily
	limited by climatic factors, recent developments of cSDMs have tried to incorporate non-climatic
	factors such land use, interspecific interactions or dispersal abilities to improve the modelling
	tools

ACKNOWLEDGEMENTS

We are grateful to Georges Kunstler, Nicolas Martin-St Paul, Cédric Perret, Estelle Barbot, Volker Bahn and Maxime Dubart for useful discussions and comments, Quentin Legros for his help in programming a first version of the diffusion model and Timothée Bonnet for his statistical advice.

We also thank two anonymous reviewers for their comments on previous versions of this manuscript. This study has been funded by the CNRS and by ANR project BioProFor (contract no. 11-PDOC-030-01).

DATA ACCESSIBILITY STATEMENT

The data supporting the results will be archived in Dryad and the data DOI will be included at 884 the end of the article. The database and script have also been submitted with the present MS in Supporting Information (Metadata S1 & S2).

SUPPORTING INFORMATION 922

Additional Supporting Information for this article, including supplementary tables and figures, may be found online.