Ahmed El-Gabbas

Ahmed ElGabbas

Nature Conservation Sector

Gebel Elba National Park

Nature & Science Foundation

BioMAP Project

Albert-Ludwigs-Universität Freiburg

University of Freiburg

Department of Biometry and Environmental System Analysis

MRes in conservation Biology

University of Nottingham

Suez Canal University

The British Chevening scholarship (2011/2012)

GERLS scholarship

German-Egyptian Research Long-term Scholarship

Deutscher Akademischer Austauschdienst

Ahmed El-Gabbas

Ahmed El-Gabbas

Ahmed El-Gabbas & Francis Gilbert
2016
The Desert Beauty Calopieris eulimene: a butterfly new to Egypt (Insecta: Lepidoptera)
Zoology in the Middle East 62:3(279-281)
10.1080/09397140.2016.1202984
The butterfly genus Calopieris Aurivillius, 1898 has only one species: the Desert Beauty Calopieris eulimene (Klug, 1829) (Ackery, Smith, & Vane-Wright, 1985), which is a rare Afrotropical species considered to be one of the most xerophilic butterflies in Africa (T. Larsen, pers. comm.; Williams, 2015). It is endemic to the Somali sub-region (Nazari et al., 2011), having been recorded only from Sudan, Yemen, western Saudi Arabia (20 km south of Mecca), Chad, Eritrea, and Ethiopia (T. Larsen, pers. comm.; Ackery et al., 1985; Gabriel, 1949; Larsen, 1982; Williams, 2015). The type specimens come from Ambukol, Dongola district, Sudan (Longstaff, 1913). The distribution is shown in Figure 1. Most of the available records are from between 1828 and 1980, with only one relatively recent record in 2007 from South of Jebel Aulla, Sudan (Williams, 2015). It is a poorly known butterfly with relatively few records and hardly any information on its biology and ecology (T. Larsen, pers. comm.; but see Waterfield, 1925).
The larvae turn up on the leafless bushes of the Desert Caper Capparis decidua (Capparaceae), and adults are mostly associated with it (Ackery, Smith, & Vane-Wright, 1985; Longstaff, 1913): it does not seem to visit the flowers of other plants (Waterfield, 1925).
Many (20-30) individuals of C. eulimene were recorded on 29 May 2011 and 29 November 2012 in a small wadi named ‘Srob Kwan’ in the Meisah area of the Gebel Elba Protected Area (22.31°N, 35.60°E; 451 m a.s.l.) in southern Egypt. The adults were found hovering over the Desert Caper (Capparis decidua; local name: Tundob; Figure 2), confirming previous observations (e.g. Ackery et al., 1985). The surrounding wadis were roughly surveyed on the same days, but no other populations of C. eulimene were observed. Other recorded butterfly species included Danaus chrysippus, Pontia glauconome, Colotis danae, C.chrysonome, and C.liagore.
Based on the two most recent comprehensive studies on the Egyptian butterflies, there are 61 butterfly species recorded from Egypt (Gilbert & Zalat, 2007; Larsen, 1990). C.eulimene has not been reported from Egypt before, although Larsen (1990) expected its distribution to extend to the extreme south of Egypt.
The Desert Caper is widespread in Egypt, inhabiting desert wadis and sandy alluvial plains. It is found in various phytogeographical regions in Egypt, including the Nile region, oases, desert areas, the Red Sea coastal strip, Sinai, and Gebel Elba (Boulos, 1999). In Gebel Elba Protected Area, it has been recorded from many locations (for details: Al-Gohary, 2009). It has also been recorded from other locations in southern Egypt, including Wadi El Gemal and Wadi El Allaqi Protected Areas (Mahmoud & Gairola, 2013; Springuel, Sheded, & Murphy, 1997). The presence of the Desert Caper at different locations in southern Egypt may support the existence of C. eulimene populations at other places in southern Egypt, especially since such places are often remote areas that receive very low sampling effort: sampling is commonly biased towards areas close to main cities or roads. Given that C. eulimene has been recorded so far only from a very small region in Egypt, it gives the Gebel Elba Protected Area high responsibility to conserve this species and its habitat, and other nearby Protected Areas (Wadi El Allaqi, Wadi El Gemal) should be thoroughly explored to show more accurately its distribution within Egypt.
Acknowledgements: We are grateful for the late Dr. Torben Larsen (†) for helping with the identification of C.eulimene. Sincere thanks for Mahmoud Saber and Ahmed Badry for assistance in the field. We also thank the African Butterfly Research Institute, Nairobi for providing location details of the Saudi Arabian record. Figure 1. Distribution of the Desert Beauty (Calopieris eulimene). Blue dots show its known distribution from Sudan, Yemen, Ethiopia, Eritrea, and western Saudi Arabia. The area shaded with gray in southeastern Egypt shows the location of the Gebel Elba Protected Area, where C. eulimene were observed.
Figure 2. Panoramic view of Wadi Srob Kwan, where the specimens of Calopieris eulimene were recorded (above) and a bush of the Desert Caper Capparis decidua, the food plant of the larvae of Calopieris eulimene.

Ahmed El-Gabbas
Sherif Baha El Din
Samy Zalat
Francis Gilbert
Conserving Egypt’s reptiles under climate change
Journal of Arid Environments, 127:211–221
10.1016/j.jaridenv.2015.12.007
Climate change has caused range shifts and extinctions of many species in the recent past. In this study, the effects of climate change on Egyptian reptiles were investigated for the first time using species distribution models. Maxent was used to model the current and future distributions of suitable habitats for 75 terrestrial reptile species from Egypt. The modelled distribution for current suitable conditions for each species was projected into the future at three time slices using two emission scenarios from four global circulation models and under two assumptions of dispersal ability. Climate change is expected to vary in its effects spatially, with some areas characterized by increased species richness while others show declines. Future range changes vary among species and different future projections, from the entire loss to large gains in range. Two species were expected to become extinct in at least one future projection, and eight species were expected to lose >80% of their current distribution. Although Protected Areas have greater conservation value, on average, compared to unprotected areas, they appear inadequate to conserve Egyptian reptiles under expected climate change.
Over the last century, relatively rapid changes of the Earth's climate have been recorded with warmer temperatures accompanied by altered geographical and seasonal distributions of precipitation (Araújo and Rahbek, 2006; Thuiller, 2007). There is widespread agreement that climate change will have a large impact on the survival of populations, species and communities (Suarez and Tsutsui, 2004), and that biodiversity is continually being transformed in response to it (Hannah et al., 2005). Over the last 40 years climate change has been implicated as the main cause for distributional shifts and extinctions (Thomas et al., 2004), with a particularly strong impact on butterflies, birds and species at high altitude (Hannah, 2011). Although the recorded effects of climate change on biodiversity seem in general to materialize rather slowly, its effects are expected to become increasingly prominent over the next 50 years and beyond (Thuiller, 2007). Some climate change model forecasts suggest that 15–37% of current species will be committed to extinction by 2050 (even without taking into consideration the biological factors of competition and evolutionary history) (Thuiller, 2007), making it essential to involve measures for mitigating its potential impacts in future conservation plans.
Detailed information on the ecological and geographical distributions of species is essential for conservation planning and forecasting (Elith et al., 2006) especially where species face multiple conservation problems. Species Distribution Models (SDMs) quantify the relationship between species' occurrence and environmental variables and allow users to extrapolate this relationship to new areas or time periods. SDMs have been widely used to estimate the potential impacts of climate change on species distributions and ecosystems (Franklin, 2009) and estimate potential future extinction risks. Once a model has been calibrated for current climate conditions, it can be used to estimate potential distributions at different time periods (in the future or the past) by using information on expected climates, or to different study areas in order to assess the potential locations where invasions are more likely to establish (Franklin, 2009). This helps to manage species facing possible future threats by identifying biological corridors for dispersal, determining sites for re-introduction and areas require more protection measures (Thuiller, 2007).
Climate change will potentially affect the biodiversity and species composition of Egyptian ecosystems (Tolba and Saab, 2009), but only a handful of quantitative studies of the local fauna/flora have been conducted using SDMs. This may be because the models are relatively new and because until very recently reliable comprehensive biodiversity records for the Egyptian fauna and flora were not available. As a developing country, Egypt lacks a recording scheme to collate species sightings (either at national or even at local Protected Area level). Between 2004 and 2008, all available biodiversity records were collated by the BioMAP project (see: http://www.biomapegypt.org/), and studies of butterflies and mammals using SDMs became possible (Gilbert and Zalat, 2008; Basuony et al., 2010; Newbold, 2009; Newbold et al., 2009a,b). A few other SDM studies have been published: El Alqamy et al. (2010), to estimate the potential distribution of the Nubian Ibex (Capra nubiana) in South Sinai; Soultan (2011), to test the potential impact of climate change on the distribution of four Egyptian antelope species; and Leach et al. (2013), who assessed the effect of climate change on the future effectiveness of the Protected Area network, using data on Egyptian butterflies and mammals. The only other study of the effects of climate change in Egypt, by Hoyle and James (2005) on the world's smallest butterfly, the Sinai Baton Blue (Pseudophilotes sinaicus), used an occupancy-based Population Viability Analysis.
The global current Protected Areas seem not to overlap well with areas of high biodiversity value, and are traditionally determined spatially and environmentally under the assumption of relatively low changes in species distribution in the future (Araújo, 2009; Leach et al., 2013). As climate change is expected to affect the future distribution/range of many species globally (potentially moving some species out from Protected Areas; Hannah et al. 2007), it challenges the future effectiveness of current Protected Areas. Future conservation investments should be effectively prioritized due to the limited resources available (Wilson et al., 2009; especially in the developing countries), and early actions may be more effective and less costly than delayed or no actions (Hannah et al. 2007). Conservation prioritization can be at taxonomic (for the protection of some rare or endangered species) or spatial (conserving a particular habitat type or species rich areas, e.g. potential Protected Areas) scale. Spatial conservation prioritization uses spatial quantitative data to identify areas of high conservation priority (Wilson et al., 2009), and some techniques have been available recently; e.g. Zonation (Moilanen et al., 2012) and Marxan (Ball et al., 2009). Some of those techniques can use spatial outputs from SDMs to prioritize areas for conservation under current and future climates.
To date, there are 30 Protected Areas in Egypt covering >15% of its total area. Their distribution shows good spatial coverage, although some areas of high diversity importance (especially at the Nile Valley and Delta) are not yet protected (Newbold et al., 2009a). They were declared relatively recently (first in 1983) and were determined mainly based on experts' known knowledge of the country biodiversity (Newbold et al., 2009a). The capacity of current Protected Areas in Egypt to mitigate for potential impacts of climate change on different species groups is not well-investigated yet and a qualitative assessment of their future effectiveness is highly required. Also, spatial prioritization of the Egypt's landscape (inside and outside the Protected Areas) is required to identify potential locations for future Protected Areas and identify current Protected Areas need more conservation effort in the future.
As a developing country, Egypt lacks enough good-quality data to be used directly for spatial conservation prioritization, but SDMs can provide valuable estimates for species suitabilities in the space. In this study, we use data for the Egyptian reptiles, for the first time, to run SDMs (as a representation group for the Egyptian fauna). Baha El Din (2001, 2006) presented a geographic interpretation of Egypt's herpetofaunal distribution, qualitatively identifying priority conservation areas, but very little has been published on how the Egyptian herpetofauna may respond to climate change. We used Maxent to model the distribution of Egyptian reptiles under current climate conditions, then models are projected into the future to show how collectively and individual species are expected to respond to future climate change under different assumptions. For each species, future expected range change (% loss or gain of currently suitable habitats) is estimated, aiming to shed light on some species require more strict conservation actions. Expected reptiles' species richness (under current and future climates) is estimated to identify areas of current high reptile suitability and areas expected to undergo much changes in suitability in the future. Model predictions were used for prioritizing the Egyptian landscape under current and future climates. We used Zonation software (Moilanen et al., 2012) to assess the likely effectiveness of Egypt's Protected Area network under current and future climate. Outputs from Zonation represent hierarchic ranking of the landscape for conservation and can be easily visualized as maps. We hope the results of the current study to be useful for biodiversity conservation in Egypt and (along with the results of relative studies) to be taken into account by the decision makers in future national conservation plans.
Study area and presence records
According to Baha El Din (2006), Egypt has at least 109 species of reptiles: 61 lizards, 39 snakes, a crocodile, seven turtles and a tortoise, and Acanthodactylus aegyptius has been added as a distinct species since (Baha El Din, 2007). The main source of species location records was the BioMAP (www.biomapegypt.org) database of Egyptian biodiversity records; with our own records, the data are comprehensive for the herpetofauna. Records for marine or aquatic species from the Nile were excluded because of the lack of GIS predictor layers for aquatic environments.
The taxonomic and georeferencing accuracy of these records were exhaustively checked, assessed and revised. We limited our scope to species with at least eight unique pixels (at a resolution of 2.5 arc minutes) to avoid over-fitting (Baldwin, 2009): this meant that a total of 75 species were processed (49 lizards, 25 snakes, and a tortoise: Table S1). The coverage of the records is good (Fig. S1a); they include most of Egypt's landscape and habitats. There was an inevitable bias in recording effort (represented by the number of valid records) towards the main cities and populated areas (Fig. S1b): unsurprisingly, the highest collection effort was found around the greater Cairo district, followed by South Sinai (the St Katherine area), the Alexandria area, some areas around Fayoum and Wadi El-Natrun and small patches near El-Arish and Mersa Matruh. (Figs. S2–S3 show a map of Egypt overlaid with geographical locations of areas mentioned in this study and the locations of currently established Protected Areas, respectively).
Environmental predictor variables
Climate data for the near past (1950–2000) were downloaded from WorldClim Global Climate Data v1.4 (release 3 – see http://www.worldclim.org/bioclim) (Hijmans et al., 2005). These data represent current climate conditions and comprise 19 ‘bioclimatic’ variables derived from precipitation and temperature records (Table 1). Four (Bio-7, -14, -17 & -18) were rejected a priori as not sufficiently variable across Egypt, reflecting the country's status as the driest country in the world (FAO, 2012). A resolution of 2.5 arc-minutes (∼5 × 5 km) best matched the level of positional uncertainty that accompanies the records (most of which derive from museum records georeferenced by site name), and because the climate data for Egypt were interpolated using relatively few weather stations largely concentrated in the Nile Valley and Delta (see Fig. 2.3 in Newbold, 2009). This relatively coarse resolution is also useful to minimise the effect of ignoring species interactions on modelled distributions (Pearson and Dawson, 2003).
Table 1. List of variables used to calculate VIF values; rows shaded with grey show variables with VIF values less than 10 and so used to run the models. Four (Bio-7, -14, -17 & -18) were rejected a priori (see text).
Altitude Altitude
NDVI_Max NDVI maximum value
NDVI_Difference Absolute difference between the highest and lowest NDVI values
Bio1 Annual mean temperature
Bio2 Mean diurnal range (mean of monthly (max temp − min temp))
Bio3 Isothermality (Bio2/Bio7) (*100)
Bio4 Temperature seasonality (standard deviation*100)
Bio5 Maximum temperature of the warmest month
Bio6 Minimum temperature of the coldest month
Bio8 Mean temperature of the wettest quarter
Bio9 Mean temperature of the driest quarter
Bio10 Mean temperature of the warmest quarter
Bio11 Mean temperature of the coldest quarter
Bio12 Annual precipitation
Bio13 Precipitation of the wettest month
Bio15 Precipitation seasonality (coefficient of variation)
Bio16 Precipitation of the wettest quarter
Bio19 Precipitation of the coldest quarter
Elevation data were obtained from the SRTM Digital Elevation Database v4.1 [available at: http://www.cgiar-csi.org/data/elevation]. The tiles that cover Egypt were downloaded at 90-m resolution and rescaled to 2.5 arc minutes. We also used an Egyptian habitat map from the BioMAP project (A. A. Hassan, unpublished data). In this map, the habitat of Egypt is classified into 11 classes (sea, littoral coastal land, cultivated land, sand dune, wadi, metamorphic rock, igneous rock, gravels, serir sand sheets, sabkhas and sedimentary rocks – see: Newbold et al., 2009a). Normalized Difference Vegetation Index (NDVI) data for seven years (from Jan 2004 to Dec 2010) were downloaded from the SPOT Vegetation website (see: http://free.vgt.vito.be/). Two layers were derived from these maps and used as predictors in the models: maximum NDVI value (indicating how much vegetation there is per pixel), and the difference between maximum and minimum NDVI value (indicating variability in vegetation per pixel).
To reduce colinearity among predictors and over-fitting, the Variance Inflation Factor (VIF) was calculated for the continuous predictors using R v2.15.0 (R Development Core Team, 2012) and used to prune predictors before use (e.g. Bombi et al., 2012) until all remaining variables yielded VIF values below 10 (Table 1): these were then used in the Maxent modelling.
Climate models used
To estimate the potential impact of climate change, current distribution models were projected into the future for three time slices (2020s, 2050s, and 2080s, each the median of a 30-year window). Future climate data (IPCC, 2007: downloaded from http://www.ccafs-climate.org/) derived from four independent Global Circulation Models (GCMs) were used to obtain consensus estimations about the fate of the herpetofauna. These GCMs were the Hadley Centre Coupled Model, version 3 (HadCM3), the Canadian second-generation coupled global climate model (CGCM2-CCCma), the Australian model (CSIRO Mk2) and a Japanese model (NIES99), all frequently used in similar studies (e.g. Araújo et al., 2006; Holt et al., 2009). Our results assume that species try to track their suitable habitat (dependent on dispersal ability) rather than accommodate or adapt to the new conditions.
For each GCM model and time slice, two emission scenarios (A2a and B2a) were used, the most commonly used scenarios in climate-change assessments (Hannah, 2011). They reflect two different assumptions about the levels of CO2 emissions, linked to assumed demographic changes and socio-economic and technological developments (Marini et al., 2009). The A2a scenario assumes a lot of climate change caused by medium to high CO2 emission rates (the “business as usual” scenario), while the B2a scenario assumes relatively little change because of mitigation measures (the “moderate” scenario) (Sauer et al., 2011; Taubmann et al., 2011). There are no data on how NDVI and habitat variables might change in Egypt in the future, and so we assume they (and altitude) remain constant. Apart from the A2a scenario for 2080, nearly all the expected climatic conditions are similar to current values (see Fig. S4).
Species distribution modelling
We used Maxent v3.3.3k (Phillips et al., 2006 – see: http://www.cs.princeton.edu/∼schapire/maxent/) to run the species distribution models. Maxent performs well in comparison with alternative algorithms and shows equal or higher predictive accuracy on cross-validation (Elith et al., 2006; Hernandez et al., 2006). Crucial to our decision not to use ensemble modelling, Maxent is very robust to small numbers of training records, outperforms other algorithms when using a relatively small number of occurrence records (Franklin, 2009), and is robust against moderate location errors (Graham et al., 2008). We used ‘permutation importance’ (see Phillips and AT&T Research, 2011) to detect which predictor had the greatest influence on the model because it does not depend on the path used by the algorithm (Songer et al., 2012).
Ten replicated runs with cross-validation and default settings (except iterations = 1000) were made for each species, by iteratively partitioning the records to use 90% for training and 10% for testing (as recommended by Phillips and AT&T Research, 2011); this gives a more stable model performance (avoiding overfitting) and minimises the effects of errors and biases. The default logistic output format was chosen, i.e. related to probability of suitable conditions, ranging from 0 to 1. Logistic output has been criticized recently (e.g. Merow et al., 2013) because many users assume that it represents the probability of occurrence; to do that, however, the calculation must assume that overall prevalence is fixed at 0.5 for all species, which may not be accurate: estimating true prevalence is very difficult (Elith et al., 2011). Both logistic and raw outputs are monotonically related and rank locations identically (Elith et al., 2011), making rank-based measures such as AUC (for model evaluation: Elith et al., 2011) and thresholding (using values from predictions of the training presence locations, e.g. 10th-percentile threshold) identical whichever is used. We used logistic output for two main reasons: a) it standardises the expected habitat suitability values across species, permitting the adding together of species to create the ‘probability richness’ maps; b) it allows Maxent automatically to convert the probabilities into expected presence/absence (or, more precisely, suitable/non-suitable) habitat. For the latter we used the 10th-percentile training presence threshold (following Pearson et al., 2007). For each species (and for each of the logistic and thresholded distributions), the ten replicate distribution maps were converted to a single consensus map, by Maxent for logistic values, and manually for thresholded values (where each pixel was a consensus presence if it had presence values in more than half of the model runs, i.e. >5).
Distributions were then projected into three future time slices (2020, 2050, and 2080), in each of four GCM models (HadCM3, CCCma, CSIRO, and NIES99) and two emission scenarios (A2a and B2a), giving a total of 24 future projections. Thresholded future distributions involved two extreme assumptions about dispersal ability: unlimited and no dispersal. For each of the 24 projections, three consensus maps were produced; one using the probability (logistic) output, and two using the thresholded distributions (assuming unlimited and no-dispersal). The maps were then averaged across the different GCMs; for thresholded distributions, this involved assigning a pixel as potentially suitable if it had a value of presence in more than two GCMs. Thus there were 18 (three time slices, two scenarios, three consensus) maps per species.
2.5. Comparisons across species
Current and future maps of estimated species richness were calculated in two different ways: first, by simply adding together the average probability of suitable conditions for each pixel across the distribution maps of all species (this assumes unlimited dispersal; e.g. Gilbert and Zalat, 2008); and second, by adding together the thresholded distribution maps (under unlimited- and no-dispersal assumptions; e.g. Trotta-Moreu and Lobo, 2010). Maps of future species richness were obtained for each of the 24 future projections and then averaged across the four GCMs. Future potential changes in species richness (but not composition) were calculated by subtracting current species richness at each pixel.
At the species level, gains and losses in future distribution were calculated using current and future consensus thresholded distributions. Then the losses across all the 24 future potential distributions across all species were added together into a species loss map, showing which areas are expected to suffer most in losing species (for each dispersal assumption). The same was carried out for gains in suitability to produce a species gain map, showing which areas are expected to gain species (only for the unlimited dispersal assumption). Species turnover, an index of dissimilarity between current and future species composition (Thuiller, 2004) was calculated for each future projection and dispersal assumption (following Thuiller et al., 2005).
In order to assess future extinction risk and to determine which species may require more protection in the future, species range changes were calculated as the percentage loss (or gain) in suitable habitats, carried out for each dispersal assumption. Each species was classified into one of the following categories: Extinct (expected loss of the entire suitable habitat - 100%), Critically Endangered (loss >80%), Endangered (loss 50–80%), Vulnerable (loss 30–50%), Least Concern (loss 30%), gain 1 (gain 30%), gain 2 (gain 30–50%), gain 3 (gain 50–80%), gain 4 (gain 80–100%), and gain 5 (gain >100) (modified from Thuiller et al., 2005) (see Tables S2–S3).
Area prioritization for conservation
Several algorithms are now available to prioritize areas for conservation and conservation planning. Using the SDM predictions, Zonation v3.1 (Moilanen et al., 2012) created a nested spatial conservation prioritization network to evaluate the effectiveness and performance of current Protected Area network in Egypt and identify new sites for conservation. Zonation prioritizes the conservation value of the landscape hierarchically based on the conservation value of sites (cells) (Moilanen et al., 2012), by iteratively removing the least valuable cells from the edges of the landscape according to set rules. The last to be removed are the most important areas (Moilanen et al., 2012). We used two removal rules: ‘basic core-area function’, to identify important (or poor) locations where a single or a few species have important occurrences; and ‘additive benefit function’, to give more weight to locations with high species richness (Moilanen et al., 2012). Fragmentation of the chosen areas is minimized by an aggregation rule, and we chose to use ‘distribution smoothing’, for which an estimate of dispersal distance is required. The dispersal ability of almost all terrestrial reptile species is very limited (Cadby et al., 2010; Edgar et al., 2010). To be conservative, in this study the dispersal distance of all reptiles was set to be equal to 1 km.
Because some species are more important to conserve than others, each species was given a weight. Weights were developed from the global and national IUCN assessments, the world distribution, and distribution patterns within Egypt. Following Gilbert and Zalat (2008) and Basuony et al. (2010), we classified Egyptian reptiles according to their world distribution and their distribution pattern within Egypt. We assessed their IUCN Red List status in Egypt according to IUCN guideline and categories (IUCN Standards and Petitions Subcommittee, 2010), although inevitably incomplete. We gave each element a score (modified from Leach et al., 2013) to indicate relative importance, and then the sum of the scores gave the relative conservation weight of each species (see Tables S1 and S4). We used these values to weight the Zonation prioritization process.
For current and future climates, Zonation was run using the maps of habitat suitabilities for all species. The number of cells removed at each iteration was set to 10 (=‘warp factor’). Zonation produces ASCII raster files with ranked pixel values scaled to be from 0 to 1: values > 0.7 were considered to be of high conservation importance. These important areas were then overlaid with the current Protected Area network in Egypt to show if the network in Egypt is adequate to conserve the species in the face of climate change, and if there are areas outside the current Protected Areas boundaries that require special protection measures.
Model performance & variable importance
In terms of mean AUC values for the 10 cross-validation runs, the performance of current-distribution models was good (see Table S1), ranging from 0.78 to 0.99 with an overall mean of 0.93 ± 0.05. Only one species (Ptyodactylus siphonorhina) had a mean AUC of less than 0.8, while in 36 species it was greater than 0.95. As almost all species had a mean AUC >0.8, all models were accepted and processed for further analyses. More widespread species tended to have lower AUC values (correlations [n = 75] with: number of suitable pixels, rs = −0.85, p 0.005; extent of occurrence, rs = −0.77, p 0.005; number of pixels with records, rs = −0.433, p 0.005).
Variables with highest mean permutation importance across all modelled species were altitude (highest for 34 species), Bio4 (temperature seasonality, 10 species), and Bio13 (precipitation of wettest month, 9 species). Those with the lowest mean permutation importance, and thus not influencing the final models very much, were the difference between maximum and minimum NDVI, Bio15 (precipitation seasonality), habitat and Bio9 (mean temperature of driest quarter); the first and last of these were never the most influential variable for any species.
Species richness
The patterns of current species richness (Fig 1) were similar using either logistic or thresholded distributions, both emphasizing the lower Nile Valley, the Suez Canal area, and the coasts of the Mediterranean and Suez Gulf (Fig. S2 for locations map). The pattern is concordant with the species richness map in Baha El Din (2006), produced using records and qualitative criteria (see Fig. S5a). Under climate change, regardless of emission scenario or dispersal assumptions, several areas are expected to have increased species richness (Fig. 2), with generally greater magnitudes under the A2a scenario, while others are expected to decline (e.g. the Suez Canal area), with greater declines under the B2a scenario by 2050 and 2080.
Fig. 1. Current estimated reptile species richness of Egypt, using the sum of: (a) Maxent logistic predictions; and (b) thresholded predicted distributions.
Fig. 2. Estimated future changes in species richness as a result of anthropogenic climate change under two IPCC scenarios (A2a, B2a) and at three times in the future (2020, 2050, 2080). Each map is the difference from the current species richness (shown in Fig. 1a), and is averaged across four GCM models. Maps made (a) using Maxent logistic output; and (b) using thresholded distributions (see Methods).
Areas with the highest expected loss of species include between greater Cairo, Ismailia and Suez, south of Suez on both sides of the Suez Gulf, and Wadi El-Natrun, much greater under the B2a scenario (Fig. 3a). If species can disperse well, then the highest gains in species are expected to include a large proportion of the northern half of the Western Desert, again greater under the B2a scenario (Fig. 3b). In terms of percentage change of the reptile assemblage, if dispersal is unlimited then the Western Desert is expected to undergo high turnover (Fig. S6a), but this is because there are so few species there currently. Under no dispersal, turnover is very limited (Fig. S6b), again with highest values again in the Western desert because of the species poor assemblages there.
Fig. 3. Expected patterns of cumulative (a) losses, and (b) gains of species as a result of climate change. The losses of each species in each pixel (under both dispersal assumptions) and gains (under unlimited dispersal only) were averaged across GCM models, accumulated across species, and then plotted.
Range and status changes
With unlimited dispersal, no species is expected to become “extinct” (100% loss of suitable habitat) under all or the average of all GCMs. There are a couple of species expected to become extinct by losing their entire area of suitable habitat in at least one of the future projections: Tarentola mindiae and Hemidactylus robustus. Using the average gain or loss of suitable habitat across the four different GCMs, eight species are expected to be classified as Critically Endangered (loss > 80%), twenty as Endangered (loss 50–80%), and twenty-two as Vulnerable (loss 30–50%) in at least one mean future projection (details, see Table S2).
When assuming no dispersal at all, none of the species are expected to lose all suitable habitats under all or the average of all GCMs, but T. mindiae and H. robustus again experience complete loss of habitat under at least one of the projections. Using the average loss of suitable habitat across the four different GCMs, ten species are expected to be classified as Critically Endangered, twenty-seven as Endangered, and thirty-three as Vulnerable in at least one mean future projection (details, see Table S3).
Area prioritization for conservation
The areas with the highest current prioritization value for conservation were similar for both cell-removal methods, mainly emphasizing the north coast, Suez Canal area, South Sinai, the Gebel Elba region, Qattara Depression, Wadi El-Natrun, and around Cairo & Fayoum, with the ‘additive benefit’ function putting greater emphasis on the Nile Delta and its fringing areas (Fig. 4). This overall pattern does not change much under future climate projections, though a few areas decline in priority and much greater area of the northern Western Desert increases. These effects are more pronounced with the ‘core-area function’ used as the cell removal rule.
Current conservation prioritizationl of areas, according to the Zonation…
Fig. 4. Current conservation prioritizationl of areas, according to the Zonation algorithm under two cell removal rules: (a) core-area function (identifies areas where one or more species have important occurrences); and (b) additive benefit function (which gives greater emphasis to areas with high species richness). Darker colours represent higher conservation value. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
The mean prioritization value in all models was higher in Protected Areas than outside them (Fig. 5), with both slightly increasing in future projections. The difference in prioritization value between inside and outside Protected Areas declined with time, especially by 2050 and 2080 under the B2a scenario.
Fig. 5. Mean prioritization value for conservation (±95% confidence limits) derived from Zonation, for Protected Areas (PAs) and outside Protected Areas (non-PAs) for current and expected future conditions, using the cell-removal rule of (a) core-area function; and (b) additive benefit function.
The results of this study show for first time the potential impacts of climate change on the distribution of Egyptian reptiles. Some reptiles are expected to lose much of their currently suitability areas (up to 80% in some cases), and some areas will lose or change a large proportion of their species. Thus there will be a need for greater conservation efforts in the future. Some of these areas are already inside Protected Areas and so hypothetically they are protected, while others are either just outside or very far from Protected Areas. Although overall Protected Areas contain better (higher priority) values compared to unprotected sites, the current Protected-Area network appears inadequate for future conservation. Our results together with others (e.g. Leach et al., 2013) will help the Egyptian Environmental Affairs Agency and other decision-making parties in Egypt to direct conservation efforts and limited conservation resources in the right direction.
Model performance & variable importance
The negative correlations found between AUC and the area occupied, and with the extent of occurrence, concur with findings of many other studies (e.g. Brotons et al., 2004; Elith et al., 2006; Hernandez et al., 2006). Rare species are usually habitat specialists, with narrow environmental tolerances, and are environmentally or geographically restricted compared to widespread species (Elith et al., 2006; Pandit et al., 2009). Widespread species are likely to be generalists occupying a wide range of habitats and climates, making it difficult to distinguish between suitable and unsuitable habitats (Franklin et al., 2009). However, as pointed out by Elith et al. (2006), this may be an artefact, an inevitable result for species with small relative occurrence areas, and perhaps a reason not to use AUC to compare model performances (Lobo et al., 2008). Most (about 60%) of the Egyptian herpetofauna are narrowly distributed, occupying less than 10% of Egypt's area (Baha El Din, 2006). There is no evidence of a strong correlation between model accuracy and the number of records used (see Elith et al., 2006; Newbold et al., 2009b), and the weak negative correlation found here is unusual since most (e.g. Stockwell and Peterson, 2002; Hernandez et al., 2006) find that models with larger sample sizes have higher accuracy (but see de Pous et al., 2011).
We found altitude to be the most effective predictor for many species. In contrast, the NDVI predictors were not useful, unlike in many other cases (e.g. Egbert et al., 2002; Anderson et al., 2006; Kgosiesele, 2010; Taheri, 2010), but perhaps typical for arid environments (e.g. Soultan's (2011) study on Egyptian antelopes) or for herpetofauna (e.g. Teodoro et al., 2013). Habitat categories likewise made only a low contribution to the reptile models.
Species richness
Species richness by itself is not adequate as an indicator of biotic change, since species composition can change independently. Our results suggest that in the future some areas with currently high expected species richness may lose many species, or undergo large changes in species composition. Some of these areas are already under a certain degree of protection, but even if they can be maintained or enhanced (e.g. Siwa oasis, Gebel Elba, and the coastal areas over the Aqaba gulf: Fig. S2), current Protected Areas are not sufficient to prevent expected loss of species. Other areas (e.g. Wadi El-Natrun, the Suez Canal area, Red Sea inland wadis, the poorly known Qattara Depression, and the areas around Gebel El-Gallala and Gebel El-Hallal) are unprotected and measures are needed to conserve them. There are similar situations elsewhere: for example, high species loss and turnover in the reptiles of a priority conservation hotspot in Europe (the Iberian Peninsula) are expected (de Pous, 2011), with great variation in expected patterns of loss and change with different projections and dispersal assumptions. Other studies have also concluded that large changes are likely in chelonians worldwide (Ihlow et al., 2012), Mediterranean mammals (Maiorano et al., 2011), US trees and vertebrates (Currie, 2001) and MesoAmerican trees (Golicher et al., 2012).
Range and status changes
Two species are expected to lose their suitable habitats entirely in at least one future projection: 1) T. mindiae is a near-endemic species recorded just from northwest Egypt and northern Cyrenaica (eastern Libya); its distribution in Egypt is restricted to Siwa oasis, the Qattara Depression and their periphery (Baha El Din, 2006). Future range changes suggest it will be Endangered (loss 50–80%) by 2020 and as Critically Endangered (loss >80%) by 2050 and 2080 (under all emission scenarios and dispersal assumptions). 2) H. robustus is localized to the East African coast from Zanzibar to southern Egypt, Arabia, east to Pakistan; its distribution in Egypt is on the Red Sea coast from El-Quseir southwards (Baha El Din, 2006). The A2a scenario suggests it will be Critically Endangered by 2020 and 2050, while under the B2a scenario it will be Endangered by 2020 and Critically Endangered by 2050 and 2080. More than half of its current distribution in Egypt is located within the Wadi El-Gemal and Gebel Elba Protected Areas, so hypothetically it is protected. However, it is found mainly on buildings and hence associated with man; it may in fact be an invasive to Egypt.
Another eight species are expected to be classified as Critically Endangered (loss >80%) in at least one averaged future projection (Acanthodactylus longipes, Cerastes vipera, Eryx jaculus, Eumeces schneiderii, Malpolon moilensis, Ptyodactylus guttatus, P. siphonorhina and Trapelus mutabilis). According to the current assessment, twenty-seven species are classified as Endangered (loss 50–80%) in at least one projection (mostly by 2080 in the A2 scenario); two of these are classified, assuming no dispersal, as Endangered in all future projections (Uromastyx aegyptia and Pristurus flavipunctatus; for details see: Tables S2–S3).
Reptiles are also likely to gain from global warming, especially in more northern areas of the globe, if they can disperse and expand their distributions, but under perhaps the more realistic assumption of no dispersal, most species were expected to lose habitat suitability by 2050 (Araújo et al., 2006). In Egypt, areas of clear future expected gain (assuming unlimited dispersal) are the area between the Siwa oasis and Qattara Depression, northwards almost to the northern coastal strip, and some sparse locations in central Sinai. There were greater gains under the B2a scenario, in the eastern part of the Western Desert and some patchy areas around the Nile (Fig. 3b). This pattern is similar to that reported by Araujo et al. (2006).
Some of the species identified here as being most at risk due to future climate change, are in fact resilient species that are ecologically very adaptable and are currently experiencing range expansion and population growth, at least in Egypt (e.g. H. robustus). This disparity is probably because modelling is based solely on the probable impacts of climate change on the spatial distribution of suitable habitats, without taking evolution or ecological adaptability and resilience into account. Species vary greatly in their response to ecological change, some adapting very well even to very rapid anthropogenically induced changes (e.g. invasive species), while others are much less adaptable and tend to go extinct first.
There are no published studies that discuss the effect of climate change on Egyptian reptiles, obviating any comparisons. Information on the reptiles of the adjacent countries of Sudan and Libya are very limited, making it impossible to assess any movements that might compensate for Egyptian losses. The patterns of expected changes in Leach et al.'s (2013) similar study on mammals and butterflies do not concur with ours. Mammal species richness was expected to decline in the Mediterranean and Red Sea areas (by 40–60%) and increase elsewhere (by 80–100%); while for butterflies, declines are expected over almost the whole of Egypt (by 40–60%) except the south, which is expected to increase (by 40–60%). In the models of European reptiles (de Pous, 2011), between eight (unlimited dispersal) and 21 (no dispersal) species were expected to become extinct in at least one future projection, with two more losing more than 80% of their currently suitable habitats. This is a much more damaging scenario than the expected responses of Egyptian reptiles in our study. There are no other comparable studies from North Africa or the Middle East.
Area prioritization for conservation
In our study, Protected Areas (Fig. S3) were more valuable in conserving reptiles than areas outside them both currently and in the future, although the difference is expected to get smaller with time (Fig. 5). Although this is encouraging, the Protected Area network appears inadequate to conserve Egyptian reptiles. New Protected Areas that have already been proposed would go some way towards reducing this deficiency in coverage, although others will be needed: in the middle to north Sinai (especially around Gebel El-Hallal), the Suez Canal area, both sides of the Suez Gulf, the northern Mediterranean coast, coastal wadis of the Red Sea between Hurghada and Mersa Alam, and the Gebel El-Gallala area. More effective protection is required in the established Protected Areas, including those in South Sinai, Gebel Elba, and Siwa oasis due to their relatively high expected future suitability loss. The Protected Areas of Siwa oasis and El-Gilf El-Kebir plus the areas of Wadi El-Natrun, the Gebel El-Hallal area and wadis of the Red Sea are expected to have particularly high turnover in species composition (Fig. S6), and hence probably need more conservation efforts. Some comparable information is available for Egypt's butterflies and mammals (Gilbert and Zalat, 2008; Basuony et al., 2010; Leach et al., 2013). Egypt's network of Protected Areas appears adequate to conserve Egyptian butterfly hotspots, except for the Mediterranean coast between Alexandria and Sallum (Gilbert and Zalat, 2008), but inadequate to conserve important mammal areas, with new Protected Areas needed to be constructed in the lower Nile Valley, along the north coast between Alexandria and Sallum, top part of the Suez Gulf, and perhaps the Qattara Depression (Basuony et al., 2010). Suggestions for conserving mammals are consistent to some extent with our results, which also suggest new Protected Areas along the Suez Gulf, Mediterranean coast from Mersa Matruh to Sallum, and the Qattara Depression. Leach et al. (2013) concur with our results in suggesting priority areas along the western Mediterranean coast, around Greater Cairo, central and north Sinai, and the Red Sea coast.
Limitations
There are some limitations of the results shown in this study. 1) Egyptian climate data in WorldClim are interpolated from rather few non-randomly placed weather stations in Egypt and adjacent countries. Although making the climate data less reliable than elsewhere, there are no alternative climate data available. 2) An ensemble approach (e.g. using BIOMOD; Thuiller et al., 2009) is often suggested for SDM modelling to avoid dependence on one particular algorithm. This is not possible without a greater number of records, since other algorithms are not as flexible or robust as Maxent. Despite all the caveats with using SDMs, and criticisms of Maxent in particular (Royle et al., 2012; but see Hastie and Fithian, 2013; Elith et al., 2006; Hernandez et al., 2006), in practice countries such as Egypt have little else to guide them in attempts at climate change adaptation. 3) There were signs of sampling bias in the available records. Although this violates the assumptions of all SDMs, it is a problem common to nearly all datasets worldwide. New techniques may allow for the correction of sampling bias; e.g. using target-group approach (Phillips et al., 2009), although the data requirements are currently impossible for Egypt and most other developing countries), 4) Using Maxent's default settings is criticized by some authors: they need to be modified dependent on the data (Warren and Seifert, 2011; Radosavljevic and Anderson, 2014; Merow et al., 2013); but we think using the default settings is enough in this study.
Conclusion
Although Egyptian datasets exist, knowledge of Egyptian species and conservation priorities are limited (with consequences explored in Leach et al., 2013); nevertheless our results and those of Leach et al. (2013) encourage us to believe that useful insights can be gained for conservation priorities by using the combination of species distribution modelling plus spatial conservation prioritization techniques (Hannah, 2010; Klorvuttimontara et al., 2011). It is crucial for the Egyptian Environmental Affairs Agency to establish an accurate and reliable database of records, assessed by experts, to help improve the basic information for future conservation assessments and analyses.
Although many SDM studies have been done on a variety of organisms in recent decades aiming to solve conservation problems, evidence of their practical utility in real-world conservation is sparse, with very few studies clearly targeting conservation decisions (Guisan et al., 2013). If implemented correctly, SDMs can play a critical role in supporting spatial conservation decision-making (Addison et al., 2013). For reptiles, there is a shortage of studies estimating their response to expected climate changes (Araújo et al., 2006). Although some suggest that reptiles (and to a greater extent, amphibians) will be strongly affected by climate change (e.g. Gibbons et al., 2000; Araújo et al., 2006; Carvalho et al., 2010), very few of them were concerned with arid environments or the Middle East. We hope that this and other similar studies that estimate the fate of reptiles under climate change help decision-making authorities to take the measures necessary to conserve them.
We thank the Italian Cooperation (Debt Swap) in Egypt for funding the BioMAP project that allowed the collation of the databases used in this work; Dr. Tim Newbold for allowing the use of the data from his field trips; Dr. Mustafa Fouda (Former Director of the Nature Conservation Sector, EEAA) for facilities offered during BioMAP project; and all the BioMAP staff. This work was supported by the award of a Chevening Scholarship Masters grant to A.El-G. We are especially grateful to two anonymous referees who greatly improved an earlier manuscript.
Table S1: A list of species used in this study (with the number of records for each species), their classification according to IUCN guidelines and criteria (global and national status), and distribution status worldwide and in Egypt, models’ mean (± sd) AUC scores of 10 folds cross-validation, and the final weighting used to run priotization analyses.
Cyrtopodion scabrum
Hemidactylus flaviviridis
Hemidactylus robustus
Hemidactylus turcicus
Pristurus flavipunctatus
Ptyodactylus guttatus
Ptyodactylus hasselquistii
Ptyodactylus siphonorhina
Stenodactylus mauritanicus
Stenodactylus petrii
Stenodactylus sthenodactylus
Tarentola annularis
Tarentola mauritanica
Tarentola mindiae
Tropiocolotes bisharicus
Tropiocolotes nattereri
Tropiocolotes steudneri
Tropiocolotes tripolitanus
Agama spinosa
Laudakia stellio
Pseudotrapelus sinaitus
Trapelus mutabilis
Trapelus pallidus
Trapelus savignii
Uromastyx aegyptia
Uromastyx ocellata
Uromastyx ornata
Chamaeleo africanus
Chamaeleo chamaeleon
Acanthodactylus aegyptius
Acanthodactylus boskianus
Acanthodactylus longipes
Acanthodactylus pardalis
Acanthodactylus scutellatus
Mesalina bahaeldini
Mesalina guttulata
Mesalina olivieri
Mesalina pasteuri
Mesalina rubropunctata
Ophisops occidentalis
Varanus griseus
Varanus niloticus
Chalcides cf. humilis
Chalcides ocellatus
Eumeces schneiderii
Scincus scincus
Sphenops sepsoides
Trachylepis quinquetaeniata
Trachylepis vittata
Leptotyphlops cairi
Leptotyphlops macrorhynchus
Eryx colubrinus
Eryx jaculus
Eirenis coronella
Lytorhynchus diadema
Macroprotodon cucullatus
Malpolon moilensis
Malpolon monspessulanus
Natrix tessellata
Platyceps florulentus
Platyceps rogersi
Platyceps saharicus
Psammophis aegyptius
Psammophis schokari
Psammophis sibilans
Spalerosophis diadema
Telescopus dhara
Naja haje
Naja nubiae
Walterinnesia aegyptia
Cerastes cerastes
Cerastes vipera
Echis coloratus
Echis pyramidum
Testudo kleinmanni
Hemidactylus flaviviridis: a very densely populated species and likely to be associated with man.
Pristurus flavipunctatus: based on its narrow habitat, strict containment in natural habitats, low densities, and high vulnerability of habitat, it does fit the VU category.
Tarentola mindiae: based on its current distribution, it has quite a large range but a rather narrow habitat. It may be reasonable to be classified at higher threat level (possibly Near Threatened “NT”).
Trapelus savignii: expected to have higher threat level.
Uromastyx ornata: suggested to have a lower threat level even though its range in Egypt is very small and its population size is smaller as well.
Chamaeleo africanus: Could be of lower threat category due to its large range in Egypt (essentially the whole Nile Valley) and quite dense populations which are expanding rapidly. It also has a very large African range; even invaded parts of Greece.
Acanthodactylus pardalis: should have a greater threat level. It is practically extinct!
Chalcides cf. humilis: more recent data indicate its very widespread distribution in Egypt, it occupies a range equal or greater than that of C. ocellatus.
Eryx jaculus: should be classified at a higher threat status than E.colubrinus or at least equal to it. It is almost extinct from Egypt.
Testudo kleinmanni: in reality it may be extinct.
Table S2: Species classification according to future species range change (percentage of suitable habitats lost or gained - assuming unlimited dispersal).
Abbreviations used: Critically Endangered “CR”: loss>80%; Endangered “EN”: loss 50-80%; Vulnerable “VU”: loss 30-50%; Least Concern “LC”: loss 30 %; Gain 1: gain 30%; Gain 2: gain 30-50%; Gain 3: gain 50-80%; Gain 4: gain 80-100%; Gain 5: gain >100%.
Table: S3: Species classification according to future species range change (percentage of suitable habitats lost - assuming no-dispersal).
Fig. S1: The distribution of (a) all Egyptian reptile records; and (b) the number of records per grid square at a scale of a ¼ of a degree.
Fig. S2: A map showing the outline of Egypt’s political boundaries overlain with the main cities and geographical locations mentioned in this study.
Ras Mohamed National Park
Zaranik Protectorate
Ahrash Protectorate
El-Omayed Protectorate
Elba National Park
Saluga and Ghazal Protectorate
St. Katherine National Park
Ashtum El-Gamil Protectorate
Lake Qarun Protectorate
Wadi El-Rayan Protectorate
Wadi Allaqi Protectorate
Wadi El-Assuti Protectorate
El Hassana Dome Protectorate
Petrified Forest Protectorate
Sannur Cave Protectorate
Nabq Protectorate
Abu Galum Protectorate
Taba Protectorate
Lake Burullus Protectorate
Nile Islands Protectorates
Wadi Degla Protectorate
Siwa
White Desert
Wadi El-Gemal/Hamata
Red Sea Northern Islands
El-Gilf El-Kebir
El-Dababya
El-Salum Gulf
El-Wahat El-Bahreya
Mount Kamel Meteor Protectorate
Fig. S3: A map showing the outline of Egypt’s political boundaries overlain with the Protected Areas.
Fig. S4: Average MESS (Multivariate Environmental Similarity Surfaces; Elith et al. 2011) maps of different global circulation models showing areas of future novel climates.
Colours ranges from blue (indicating similar future climate conditions compared to the current; the darker the blue, the higher the similarity) to red (indicating dissimilar climates compared to the current; the darker the red colour, the higher the dissimilarity). Results for dark red areas should be interpreted with caution.
Fig. S5: (a) The number of recorded/suspected amphibian and non-marine reptile species per 0.5º grid (from Baha El Din, 2006); (b) the predicted number of species under current conditions from this study (from thresholded distributions).
Fig. S6a: Future species turnover (a measure of dissimilarity between current and future species composition) assuming unlimited dispersal. Colour range from grey (low species turnover – small species composition change in the future) to dark red (high species turnover – high species composition change in the future).
Fig. S6b: Future species turnover (a measure of dissimilarity between current and future species composition) assuming no-dispersal. Colour range from grey (low species turnover – small species composition change in the future) to dark red (high species turnover – high species composition change in the future).

Tim Newbold, Tom Reader, Ahmed El-Gabbas, Wiebke Berg, Wael M. Shohdi, Samy Zalat, Sherif Baha El Din & Francis Gilbert 2010
Testing the accuracy of species distribution models using species records from a new field survey
Oikos, 119(8): 1326–1334
10.1111/j.1600-0706.2009.18295.x
Species distribution models are a very popular tool in ecology and biogeography and have great potential to help direct conservation efforts. Models are traditionally tested by using half the original species records to build the model and half to evaluate it. However, this can lead to overly optimistic estimates of model accuracy, particularly when there are systematic biases in the data. It is better to evaluate models using independent data. This study used independent species records from a new to survey to provide a more rigorous evaluation of distribution-model accuracy. Distribution models were built for reptile, amphibian, butterfly and mammal species. The accuracy of these models was evaluated using the traditional approach of partitioning the original species records into model-building and model-evaluating datasets, and using independent records collected during a new field survey of 21 previously unvisited sites in diverse habitat types. We tested whether variation in distribution-model accuracy among species could be explained by species detectability, range size, number of records used to build the models, and body size. Estimates of accuracy derived using the new species records correlated positively with estimates generated using the traditional data-partitioning approach, but were on average 22% lower. Model accuracy was negatively related to range size and number of records used to build the models, and positively related to the body size of butterflies. There was no clear relationship between species detectability and model accuracy. The field data generally validated the species distribution models. However, there was considerable variation in model accuracy among species, some of which could be explained by the characteristics of species.
Species distribution models have great potential as tools for conservation. They relate a set of records of the occurrence of a species to a set of variables describing relevant aspects of the environment in order to predict its distribution over the whole of the study area in question (reviewed by Wintle et al. 2005).
There is a vast amount of data on the distributions of species in museums, natural history collections and the literature (Graham et al. 2004). However, there are several limitations associated with data from these sources. First, records are often accompanied by an imprecise description of the locality from which they were taken. This translates into poor locational accuracy when the record is georeferenced (i.e. when it is assigned geographical coordinates, Graham et al. 2004). The effect of locational error in the species records on the accuracy of distributional models is generally low, but varies among different types of model (Graham et al. 2008).
Second, museum data are often biased. Such bias could be: 1) spatial – towards areas to which it is easy for scientists to gain access, or towards areas that are biologically interesting; 2) temporal – towards time periods when collecting was more frequent; or 3) taxonomic – towards species that are easy to detect or that are of more interest to the collectors (Soberón 1999, Hijmans et al. 2000, Reddy and Dávalos 2003).
The third major problem with data from museums and literature sources is that there are rarely data documenting places where the species is known not to exist (absence records) (Graham et al. 2004). There are modelling techniques designed to be used with datasets that consist only of presence records, such as climate envelope approaches and techniques that model the presences with reference to the background environmental conditions (Wintle et al. 2005). However, several of the most popular modelling approaches, such as generalized linear models (GLMs) and generalized additive models (GAMs), can only be used with both presence and absence data. A commonly-used solution to this problem is to generate random ‘pseudo-absence’ data in grid cells without presence records (Ferrier and Watson 1997). One obvious problem with using pseudo–absence data is that some absences are likely to be found in areas that are suitable for, and even inhabited by, the species (Graham et al. 2004). Of course, recorded absence of species may also prove to be erroneous. Many species are very difficult to detect and it may take many visits to a site before species absence can be inferred with any degree of confidence (Kéry 2002, MacKenzie et al. 2002). Given accurate species records from a well-designed survey, models built with only presence records have been shown to perform as well as models built with both presences and absences (Wintle et al. 2005) and may present the safest option when there is uncertainty over the reliability of absence data.
Data from museums, collections and the literature are too valuable a source of data to ignore. However, given the potential biases and inaccuracies associated with them, it is particularly important to test the accuracy of model predictions. The simplest way to test the accuracy of a species distribution model is to test its ability to predict correctly the data used to build it in the first place (Fielding and Bell 1997). This is effectively a measure of goodness-of-fit of the model. The main drawback of this approach is that a model can fit the data used to build it very well without capturing the species’ real response to the environmental variables (a phenomenon known as overfitting), and this method of model evaluation tends to lead to over-optimistic measures of model accuracy (Chatfield 1995). A better approach is to partition the data in some way, building the model with part of the dataset and evaluating it against the remainder (Fielding and Bell 1997). This is the approach taken by most studies (Hernandez et al. 2008, Franklin et al. 2009). A problem with data-partitioning approaches is that if the same bias in the species data is present in all partitions, then the model may be biased and the estimate of model accuracy inflated (Chatfield 1995). Ideally models should be evaluated using new, independent data on species occurrence (Chatfield 1995). With the wide availability of global positioning systems (GPS), records can be assigned geographical coordinates on collection, eliminating the problem of locational errors. Bias should be reduced as much as possible, particularly bias in environmental space (Wintle et al. 2005). Few studies have used independent data to validate models because collecting such data can be impractical, time-consuming and costly (Wintle et al. 2005, but see Loyn et al. 2001, Elith 2002, Ferrier et al. 2002, Elith et al. 2006, Williams et al. 2009). To the best of our knowledge only one study, on Mexican birds (Feria and Peterson 2002), has used new, independent records to test the accuracy of distribution models based on museum data. Given the potential limitations with records from museums, it is particularly important that the accuracy of models based on them are evaluated rigorously. We used independent records to test the accuracy of distribution models for a variety of species in three separate taxonomic groups.
Even if one is confident of a lack of bias in the data, different kinds of species may be more or less suited to the model-building process. There have been attempts to assess differences among species in the accuracy of their distribution models (Kadmon et al. 2003, Berg et al. 2004, Seoane et al. 2005, Hernandez et al. 2006, Newbold et al. 2009b). These studies have often found that species that are more narrowly distributed produce more accurate distribution models, possibly because small-ranged species have better-defined habitat requirements and tend to inhabit a greater proportion of the suitable environment, or because in species with larger ranges, populations show local adaptation to the environment in different areas (Stockwell and Peterson 2002, Brotons et al. 2004, Segurado and Araújo 2004, Hernandez et al. 2006, Newbold et al. 2009b). On the other hand, effects of range size could be a statistical artefact associated with the use of pseudo-absence data (Lobo et al. 2008).
Species that are easier to detect are likely to have more complete occurrence data. This may result in more accurate distribution models for these species (Seoane et al. 2005). For example, Pöyry et al. (2008) showed that the accuracy of distribution models for butterflies was positively related to wingspan, possibly owing to differences in detectability during surveys.
In this study, we modelled the distributions of Egyptian butterfly, mammal, reptile and amphibian species using records from museums, collections and the literature, presenting a rare test of their accuracy using new, independently-collected survey data as well as a test using the more traditional data-partitioning method. It was not possible to collect new species records systematically or randomly in the time available because of the remoteness and inaccessibility of many parts of Egypt, but the records were completely independent of the data used to build the models, were designed to be representative of as many habitat types as possible, given the constraints imposed by the logistics of sampling in a remote and hostile environment, and were georeferenced using a GPS and so had negligible locational error. We used the new survey data, which contain both presence and absence records, to test whether a negative effect of species range size on model accuracy persists in the absence of statistical artefacts. We also tested whether model accuracy is related to species detectability and body size (of butterflies).
Data and methods
Distribution models were compiled for Egyptian butterfly, mammal, reptile and amphibian species using Maxent ver. 3.1.0 (Phillips et al. 2006). Maxent uses a machine-learning process to produce a model where the frequency distribution of modelled probabilities is as close to uniform as possible, with the constraint that the expected value of each environmental variable (the sum, across all grid cells, of the product of the probability of occurrence and the value of the environmental variable) must equal the mean value at the presence points (the empirical average). To prevent overfitting a process called regularization is adopted, relaxing this constraint such that the expected value of each environmental variable may fall within a defined margin around the empirical average (Dudík et al. 2004). Maxent is particularly suited to use with museum data, because it designed to deal with datasets consisting only of presence records. The environmental conditions in a sample of cells from throughout the whole study area is used for comparison with the environmental conditions in cells with species presence records in.
The species data used to build the models (hereafter referred to as the original species records) were taken from the database complied as part of Egypt's BioMAP project ( www.biomapegypt.org for more details). Records were taken from museum and personal collections, and from the literature (Osborn and Helmy 1980, Larsen 1990). The identification of all extant specimens was checked by experts (Samy Zalat, Sherif Baha-El-Din, Francis Gilbert, Dr Mohammad Basuony, Al Azhar Univ., Cairo), according to the latest taxonomic opinion (Larsen 1990, Wilson and Reeder 2005, Baha El Din 2006). All records were mapped as accurately as possible using a gazetteer developed by the BioMAP project. We estimated the maximum error associated with each sampling location using the point radius method (Wieczorek et al. 2004) and removed records from highly inaccurate localities. Given positive spatial autocorrelation in the environmental variables, a moderate degree of inaccuracy in the location of species records probably does not have a large effect on model accuracy (Graham et al. 2008). The number of records available for each species ranged from 10 to 412 (median = 58); most records were made in the second half of the 20th century (Newbold et al. 2009a).
The environmental variables used in the models consisted of climate and habitat variables. The climate variables were taken from the WorldClim ver. 1.4 dataset (Hijmans et al. 2005). This dataset includes 20 variables describing altitude, temperature and precipitation (Supplementary material Appendix 1). The habitat variable used was a geological habitat classification with 11 categories (sea, littoral coastal land, cultivated land, sand dune, wadi, metamorphic rock, igneous rock, gravel, serir sand sheet, sabkha and sedimentary rock). This map was compiled using satellite imagery, and was verified by extensive ground-truthing (A. Hassan unpubl.). Using 20 environmental variables might cause the models to be overfitted (Chatfield 1995). However, the Maxent model uses a process called regularization to reduce the chance of overfitting and previous studies have shown that it can produce accurate models with small numbers of species records and similar numbers of environmental variables to our study (Wisz et al. 2008). To test whether this was the case in our study, we built a separate set of distribution models using habitat and three principal component axes based on the climatic variables. Models developed using the two sets of variables were very similar (at 5000 random points in Egypt modelled probabilities of occurrence correlated positively – Spearman's rank correlation: mean rs= 0.798 ± 0.02) and were not significantly different in accuracy (Wilcoxon's matched-pairs test: p = 0.122).
To create a second set of species data (hereafter referred to as the independent species records) with which to evaluate the distribution models, we conducted a survey of butterflies, mammals, reptiles and amphibians in Egypt in the summers (May–July) of 2007 and 2008. The reptile, amphibian and mammal species surveyed are active throughout the summer months. The flight periods of all of the butterfly species surveyed encompassed the whole period of sampling. The new records were not used to build distribution models, only to evaluate them. The new data were biased towards roads. The terrain in Egypt makes it almost impossible to sample completely randomly, with many areas situated hundreds of kilometres from the nearest road. We minimized bias in environmental space as much as possible by selecting sites that covered: (1) as large a geographical area as possible; and (2) as many different habitat types as possible, defined using a geological habitat map (A. Hassan unpubl.) and a vegetation land cover map, derived using data from the Advanced Very High Resolution Radiometer (Hansen et al. 2000). At each site we performed four 1-km walking transects at different times of day (early morning, late morning, late afternoon, evening), paced to take approximately an hour and a half each. At the same time, some members of the expedition actively searched for species in the area surrounding the start point of the transect. Transects were located such that they sampled all the major habitat types present at each site. A species was recorded as being present if it was observed at least once, and absent otherwise. Twenty-one sites were surveyed in this way (Fig. 1, Supplementary material Appendix 1). In addition to records from the fully-surveyed sites, we also included incidental observations of species from 13 other localities (Fig. 1, Supplementary material Appendix 1). Data from these sites consisted of records of species presence only, because we did not carry out replicate transects at these sites and thus could not infer species absence. Almost all new sites were situated at least 1 km from sites with records in the original dataset (Fig. 1). All fully-surveyed sites were at least three km from the nearest other site, and all but four were at least ten km from the nearest other site. Including locations with incidental records, distances among sites were sometimes much smaller; four sites were less than one km from the nearest other site and 15 sites were less than ten km from the nearest site. Butterflies were sampled by visual searching and sweep netting, reptiles and amphibians by visual scans and active searches, and mammals mainly by checking for tracks and signs, although sightings of species were also noted. Sixty species were recorded in total, 34 of which were recorded at least twice: 20 reptiles and amphibians, ten butterflies and four mammals (Supplementary material Appendix 1).
Sites with reptile, amphibian, butterfly and mammal records in the BioMAP database (grey crosses and asterisks), and sites that were sampled during the new survey (black triangles).
Imperfect detectability of species is likely to have an impact on the reliability of data describing species absence from surveys such as ours (Kéry 2002, MacKenzie et al. 2002). We modelled the detectability of species in our new survey data, following MacKenzie et al. (2002). The four transects undertaken at each site were treated as independent visits (n1, n2, n3, and n4). The likelihood (L) of obtaining a particular pattern of occurrence for a species across all four transects at all fully-surveyed sites is:
where ψ is the probability that a species occurs at a given site, p is the probability that the species is detected during one transect given that it occurs at the site, t is the transect number, n. is the number of sites where the species was recorded in at least one transect, and N is the total number of sites visited (MacKenzie et al. 2002). The parameters p and ψ were estimated using a maximum likelihood approach with the package ‘mle’ in R (R Development Core Team 2004). Upper and lower bounds of 0.0001 and 0.9999 respectively were set for both parameters. The model has been shown to be reasonably robust to sample sizes as small as those encountered here (Wintle et al. 2004). The model assumes that occurrence and detection probabilities are constant across sites, which is almost certainly not true. The modelled probabilities should therefore be considered rough estimates to gauge the reliability of the occurrence data and not as accurate estimates of the probabilities of detection and occurrence.
The distribution models were evaluated using three different sets of data. First, using partitioned data, whereby the original species records were divided randomly before modelling – half for model building and half for model evaluation. Models were evaluated using these reserved presence records and 2500 pseudo-absences (Ferrier and Watson 1997), drawn randomly from cells that lacked a record of the species in question. Second, using the independent species presence records and 2500 pseudo-absences, generated as before. Third, using the independent presence and absence records. Model accuracy was measured using the AUC statistic (Fielding and Bell 1997). The receiver operating characteristic curve is a plot of the proportion of presence records correctly predicted present (sensitivity) against the proportion of absence records incorrectly predicted present (commission) for a range of thresholds used to divide predicted presence from predicted absence. The area under the curve (AUC) measures the ability of the model to discriminate recorded presences from recorded absences (Fielding and Bell 1997). An AUC value of 1 indicates perfect discrimination and an AUC value of 0.5 indicates a model that is no better than random. Estimated accuracy according to AUC values was compared among the three approaches. We correlated estimates of accuracy made by partitioning the original species records with estimates made using the independent records, to test whether models were ranked similarly. To provide an alternative measure of accuracy to the AUC statistic, the models were also tested against the independent presence and absence records using the slope of the relationship between model predicted probability and species occurrence (presence or absence), fitted using a generalized linear model with binomial errors (McCullagh and Nelder 1989).
We tested a number of factors that may explain variation in model accuracy (measured using the independent presence and absence records) among species: (1) estimated species detectability (2) range size in Egypt; (3) number of presence records used to build the models; and (4) taxonomic group (mammals, butterflies, or reptiles and amphibians). The proportion of Egypt's land area predicted by the distribution models to be occupied was used as an index of range size. To calculate this, we converted the continuous prediction of probability of occurrence into a binary prediction of presence or absence, by assigning a threshold probability of occurrence to the model for each species. The threshold was set such that 95% of the presence records used to build the models were predicted correctly as being present (Pearson et al. 2004).
The effect of estimated species detectability on distribution-model accuracy was tested by a simple correlation test, because detectability could not be estimated for species that were not recorded during the walking transects. As an additional test of the effect of estimated species detectability, we also correlated butterfly wingspans (wing-tip to wing-tip; Gilbert and Zalat 2007) with model accuracy. The remaining factors were tested using generalized linear models with normal errors. AUC values were entered as the dependent variable, taxonomic group as a factor, and predicted range size and number of presence records used to build the model as covariates. We used a model selection method based on the approach recommended by Burnham and Anderson (2002). We built a global model with all terms, and candidate models with every combination of terms. AIC scores were extracted for each model and the difference between a model's AIC value and the lowest value of all models (the AIC difference, Δi) was calculated. Model weight was calculated using the following formula (Burnham and Anderson 2002):
where Δi is the AIC difference of the model in question and Δrs are the AIC differences of the other models. The relative importance of each variable was assessed by summing the AIC weights of all candidate models containing it (Burnham and Anderson 2002), hereafter referred to as ‘sum of AIC weights’.
Results
Estimates of the probability of detecting a species in a single transect (p) ranged from less than 0.001 to approximately 0.75 (Table 1). For butterflies, the migratory species Vanessa atalanta and Vanessa cardui, and the skipper Pelopidas thrax had low probabilities of detection, but most species were relatively easily detected. Mammals generally had much lower probabilities of detection than butterflies; the gazelle Gazella dorcas was an exception because its presence could be reliably ascertained by tracks and faeces. Reptiles and amphibians were highly variable in their estimated detectability. The snakes and the chamaeleon Chamaeleo africanus had very low probabilities of detection, while the lizards, skinks and amphibians generally had higher probabilities. Estimates of the probability of site occupancy (ψ), which is equivalent to the proportion of sites predicted to be occupied, were consistent with estimates of range size derived from the species distribution models (Spearman's rank correlation test: rs= 0.453, n = 23, p = 0.03).
Table 1. Estimated probabilities of occurrence (ψ) and detection, given occurrence (p) for species recorded in the walking transects at the fully-surveyed sites. Each transect was treated as an independent sampling event. ψ and p were estimated using a maximum likelihood approach (MacKenzie et al. 2002), assuming that both probabilities are constant across sites.
Acanthodactylus boskianus
Acanthodactylus scutellatus Cerastes cerastes
Chamaeleo africanus
Malpolon monspessulanus
Mesalina guttulata
Natrix tessellata
Ptychadena mascareniensis
Rana bedriagae
Sphenops sepsoides
Trapelus mutabilis
Butterflies
Colias croceus
Danaus chrysippus
Lampides boeticus
Leptotes pirithous
Pelopidas thrax
Pieris rapae
Pontia glauconome
Vanessa atalanta
Vanessa cardui
Zizeeria karsandra
Mammals
Capra nubiana
Cazella dorcas
Lepus capensis Model accuracy estimates made by partitioning the original species records into model-building and model-evaluation datasets, were high and significantly better than random (one sample t-test: t = 22.0, DF = 33, p 0.001). AUC values ranged from 0.666 to 0.975, with an average of 0.845 ± 0.016. Accuracy estimates made using the independent presence records (i.e. records from the new survey) and pseudo-absences were also high and significantly better than random (t = 16.7, DF = 33, p 0.001). AUC values ranged from 0.485 to 1.000, with an average of 0.875 ± 0.022. Finally, accuracy estimates generated using the independent presences and absences were reasonably high and significantly better than random (t = 4.03, DF= 33, p 0.001), although lower than estimates made using pseudo-absences. AUC values ranged from 0.219 to 1.000, with an average of 0.655 ± 0.039 (for examples of the distribution models, Fig. 2). Testing the accuracy of models against the independent records, using the slope of the relationship between model predicted probability of occurrence and observed occurrence (presence or absence), also showed the predictions to be reasonably accurate. The relationships were positive for 26/34 species, although only nine were significantly positive (GLM: p 0.05). Slope coefficients ranged from –5.67 to 22.13; the average coefficient was significantly greater than zero (one sample t-test: t = 3.16, DF = 32, p = 0.003). Estimates of accuracy made using subsets of the original presence records correlated significantly and positively with estimates made using the independent records (Spearman's rank correlation: rs= 0.544, n = 34, p = 0.001; Fig. 3).
Predicted distributions and independent occurrence records for two species: (a) the Montpellier snake Malpolon monspessulanus, which had the most accurate distribution model; and (b) the cape hare Lepus capensis, which had the least accurate distribution model. Distribution models were built with Maxent ver. 3.1.0 using records from the BioMAP database and variables describing climate and habitat. Light shading indicates areas with a high probability of occurrence, while dark shading indicates a low probability of occurrence. The independent occurrence records (+= presence; O = absence) were collected during a new field survey of 21 sites in the summers (May–July) of 2007 and 2008; these records were used to evaluate the distribution models.
The relationship, for 34 species of Egyptian mammal, butterfly, reptile and amphibian species, between distribution-model accuracy estimated using independent species presence and absence records, and distribution-model accuracy estimated using partitioned data, whereby the original species presence records were randomly divided in half for model building and model evaluation respectively. Accuracy was estimated using the AUC statistic (Fielding and Bell 1997).
Model accuracy showed no clear relationship with estimated species detectability (Spearman's rank correlation: rs=–0.294, n = 25, p = 0.154). However, for butterfly species, wingspan correlated positively with model accuracy (Pearson s correlation coefficient: r = 0.652 n = 10, p = 0.041; Fig. 4). Model accuracy was negatively related to the predicted range size of species within Egypt (GLM: sum of AIC weights = 0.952; Table 2, Fig. 5a). Surprisingly, there was also a strong negative effect of the number of species presence records used to build the models on the accuracy of predictions (sum of AIC weights = 0.991; Table 2, Fig. 5b). There was little support for an effect of taxonomic group on the accuracy of distribution models (sum of AIC weights = 0.172; Table 2).
Table 2. Results from a set of generalized linear models with a normal error distribution testing factors affecting variation in the accuracy of species distribution models among species. Factors tested were predicted range size in Egypt (R), number of presence records used to build the models (S), and taxonomic group (T). Candidate models were built with every combination of terms. These models were compared using AIC and the difference between the AIC of a model and the minimum AIC of all models. Model weights were calculated following Burnham and Anderson (2002).
For 34 species of Egyptian reptiles, amphibians, butter-flies and mammals: (a) the relationship between range size, estimated as the proportion of grid cells in Egypt predicted occupied, and the accuracy of distribution models estimated using independent species records from a new field survey; (b) the relationship between the number of presence records used to build the distribution model and model accuracy, estimated using independent species records. Model accuracy was measured using the AUC statistic (Fielding and Bell 1997).
Discussion
Overall, the distribution models built in this study were shown to be significantly better than random when tested against independent data collected by surveying a diverse range of habitats in Egypt. This strongly suggests that data from museums, natural history collections and literature can be used to make useful predictions about species’ ranges. Several studies have reached a similar conclusion (Peterson et al. 2002, Raxworthy et al. 2003), but it is rare that models are tested against independent evaluation data (but see Loyn et al. 2001, Elith 2002, Ferrier et al. 2002, Elith et al. 2006). Uncertainties and biases will be more prevalent in models built using museum and literature records (Graham et al. 2004), making evaluation with independent data more important. Some authors have experimented with using species records from separate geographical areas (Peterson and Shaw 2003, Randin et al. 2006, Heikkinen et al. 2007) or time periods (Raxworthy et al. 2003) to evaluate models. However, predictions extrapolated outside the environmental conditions encompassed by the data that were used to build the model are likely to be inaccurate in the new areas even if they are accurate in the area for which they were built. The best approach is to collect new, independent data inside the study area for which the models were developed reducing bias as much as possible, particularly bias in environmental space (Wintle et al. 2005).
The reliability of data on species absence probably depends on the relative detectability of the taxa in question (MacKenzie et al. 2002). There was substantial variation in estimated detection probability among species in the new survey. The results of the maximum likelihood model were consistent with our expectations. First, the predicted proportion of sites occupied correlated positively with predicted range size according to the distribution models. Second, detection probabilities were very low for elusive species, such as the snakes, and for rare migrants, such as the red admiral butterfly Vanessa atalanta, and higher for conspicuous and more abundant species, including some of the lizards and most of the butterflies. The accuracy of species distribution models did not appear to be affected by detection probability suggesting that, even in small-scale surveys with relatively few visits to each site, imperfect detection of species may not be a major problem. On the other hand, the accuracy of distribution models for butterfly species was positively correlated with body size, which was used as a surrogate for detectability. It is possible that our maximum likelihood-based estimates of detection probability were inaccurate; for instance, one of the major assumptions of the maximum likelihood model that we used is that occurrence and detection probabilities are constant across sites (MacKenzie et al. 2002), which is very unlikely to be true. However, very abundant and easily detectable species, such as the long-tailed blue butterfly Lampides boeticus and Bosk's lizard Acanthodactylus boskianus, had high estimated detection probabilities and inaccurate distribution models, whereas species that are difficult to detect, such as Montpellier's snake Malpolon monspessulanus, had low estimated detectability but very accurate distribution models. An alternative explanation for the relationship between butterfly wingspan and distribution-model accuracy is that larger butterflies are more mobile and able to reach a greater proportion of suitable habitat, giving a closer correlation between environmental variables and occurrence (Pöyry et al. 2008), although the effect of body size on butterfly mobility is contentious (Cowley et al. 2001).
Estimates of model accuracy made using the data partitioning approach were relatively consistent with estimates made using the new survey data. This suggests that a data-partitioning approach can give us a good idea about the relative accuracy of models and can be used to compare model accuracy among species. Accuracy estimates made using the partitioned species records and pseudo-absences, and also with independent presence records and pseudo-absences, were much higher than estimates made using both independent presence and independent absence records. This is consistent with a previous suggestion that overly-optimistic estimates of model accuracy can be generated using pseudo-absence data (Lobo et al. 2008), but it should be borne in mind that the small numbers of independent records may partly explain the low measures of accuracy using independent data. Nevertheless, further comparisons of model accuracy using pseudo-absences and real absences are needed and it would be prudent not to use data partitioning as the sole method for evaluating distribution models, especially if the models will be used for conservation decision-making. The accepted threshold of 0.7 above which models are considered to be good (Pearce and Ferrier 2000) may place undeserved confidence in poor predictions.
Some of the variation in model accuracy was explained by range size. Species with larger ranges within Egypt were modelled less accurately than species with smaller ranges. A negative effect of range size on the accuracy of species distribution models has been reported before (Stockwell and Peterson 2002, Brotons et al. 2004, Segurado and Araújo 2004, Hernandez et al. 2006, Newbold et al. 2009b), but most of these studies have used real presence data with pseudo-absence data. Thus, the apparent effect of range size could be a statistical artefact owing to pseudo-absences being more distant in environmental space from the presence records for species with smaller range sizes (Lobo et al. 2008). Our results show that the distributions of species with smaller ranges are modelled more accurately even in the absence of statistical artefacts. This could be because narrowly-distributed species have more specific climate and habitat requirements than more widespread species (Brotons et al. 2004, Hernandez et al. 2006). Alternatively, separate populations of widespread species may show local adaptation to the environmental conditions in different parts of the study area (Stockwell and Peterson 2002, Brotons et al. 2004): although two of the butterfly species have more than one sub-species in Egypt (Carcharodus stauderi and Spialia doris; Gilbert and Zalat 2007), these distinctions were not considered in this study.
Surprisingly, we found a significant negative effect of the number of species records used to build models on the accuracy of model predictions. Most previous studies have found the relationship between sample size and model accuracy, if present, to be positive (Pearce and Ferrier 2000, Phillips et al. 2004). Several studies have shown that species with narrower distributions and more specific habitat requirements are modelled more accurately (Kadmon et al. 2003, Hernandez et al. 2006, Newbold et al. 2009b). It is probable that some aspect of this was captured by sample size but not by the measure of range size that we used. For example, more narrowly distributed species are likely to be less abundant (Gaston et al. 2000) and thus detected less often during surveys. Alternatively, habitat specialists may be easier to model because they have very specific requirements, but may be restricted to particular microhabitats or resources and thus have been detected less frequently in the past.
Ideally data used to evaluate the accuracy of distribution models should be completely independent of the data used to build the models and free from any bias (Chatfield 1995), but given limited resources this may not be possible (Wintle et al. 2005). Although our new species records contained some bias (for example, towards locations near roads), we reduced environmental bias by selecting sites that covered as broad a range of climate and habitat types as possible. This approach is better than simple data-partitioning, because some bias has been eliminated and because locational error in the records has been eliminated. Moreover, it is more practicable than a truly random survey, especially for less accessible areas such as Egypt. The time constraints imposed on our field expedition meant that we were only able to survey the northeast part of Egypt. Therefore, we cannot be certain that the models were accurate for other parts of Egypt, particularly in the Western desert where museum records are scarce. Nevertheless, at least in northeast Egypt, the models appeared to provide an accurate reflection of the distribution of species.
In conclusion, our results support the use of species distribution models in ecology. Predictions made for many species in three very different taxonomic groups were shown to be accurate using completely independent species occurrence data. However, there was considerable variation across species in the accuracy of distribution models. Distribution models have great potential as tools for conservation, but it is crucial that their predictions are first evaluated thoroughly. Currently, using completely independent data to evaluate model predictions is a rare practice, which is not surprising given that conducting new surveys can be time-consuming and very expensive (Wintle et al. 2005). However, we show that even small field surveys can be used to test model accuracy and can highlight patterns in the accuracy of models.
Acknowledgements
We thank Italian Cooperation (Debt Swap) for funding the BioMAP Project; Mustafa Fouda (Director of the Nature Conservation Sector, EEAA) for facilities and comments on the work; all the BioMAP staff (Ahmed Yakoub, Alaa Awad, Muhammed Sherif, Shama Omran, Shaimaa Esa, Yasmin Safwat, Nahla Ahmed, Esraa Saber); Mohammad Basuony and Abd El Aal Attia for help during dataset preparation and preliminary analysis; and Rashed Refaey, Ahmed Refaey and other Egyptian field guides for making the field expedition possible. This work was supported by the Natural Environment Research Council (grant number NER/S/A/2006/14170).

Tim Newbold, Tom Reader, Samy Zalat, Ahmed El-Gabbas & Francis Gilbert 2009
Effect of characteristics of butterfly species on the accuracy of distribution models in an arid environment
Biodiversity & Conservation 18: 3629-41.
10.1007/s10531-009-9668-5
Species distribution models show great promise as tools for conservation ecology. However, their accuracy has been shown to vary widely among taxa. There is some evidence that this variation is partly owing to ecological differences among species, which make them more or less easy to model. Here we test the effect of five characteristics of Egyptian butterfly species on the accuracy of distribution models, the first such comparison for butterflies in an arid environment. Unlike most previous studies, we perform independent contrasts to control for species relatedness. We show that range size, both globally and locally, has a negative effect on model accuracy. The results shed light on causes of variation in distribution model accuracy among species, and hence have relevance to practitioners using species distribution models in conservation decision making.
Species distribution models have great potential to aid conservation efforts, allowing ecologists to predict species distributions over a large area, using incomplete data on species occurrence combined with maps of environmental variables, describing climate, habitat and topography (Wintle et al. 2005). Many studies have compared the accuracy of predictions made by different modelling techniques, often finding that many techniques perform similarly well (Elith et al. 2006; Hernandez et al. 2006; Phillips et al. 2006). In fact, there may be more variation in model accuracy among species than among modelling techniques (Berg et al. 2004; Elith et al. 2006). As a result, whether the characteristics of species affect the accuracy of distribution models is a question that is receiving increasing attention in the literature.
The breadth of a species’s niche has often been considered when trying to explain differences in model accuracy among species. Studies have found that species with narrow, well-defined niches are better modelled than those with broader niches (Boone and Krohn 1999; Pearce et al. 2001; Kadmon et al. 2003; Berg et al. 2004) and that models for specialist species are generally more accurate than models for generalists (Hepinstall et al. 2002; Segurado and Araújo 2004; Elith et al. 2006). Species with narrow niches generally have better-defined climate and habitat requirements, which are easier to model (Kadmon et al. 2003). Other studies have shown that the breadth of a species’ niche relative to the environmental conditions found in the study area as a whole influences model accuracy more than the breadth of a species’ niche per se (Seoane et al. 2005; Hernandez et al. 2006). More marginal species (those that have niches furthest from the average conditions of the study area) have also been shown to be modelled more accurately than less marginal species, probably for similar reasons (Luoto et al. 2005; Seoane et al. 2005; Carrascal et al. 2006; Hernandez et al. 2006). We might expect therefore that the accuracy of species distribution models will decrease with increasing niche breadth or habitat tolerance.
Models for species that have narrow distributions in geographical space have also been found to be more accurate than models for species with larger distributions (Stockwell and Peterson 2002; Brotons et al. 2004; Segurado and Araújo 2004; Hernandez et al. 2006). This may be related to the effect of niche width, with smaller range size being associated with better-defined habitat requirements (Brotons et al. 2004; Hernandez et al. 2006). Alternatively, populations of species with larger ranges may show local adaptation to different environmental conditions, decreasing the accuracy of models that consider all populations together (Stockwell and Peterson 2002; Brotons et al. 2004). Similarly, McPherson and Jetz (2007) found that endemic species were modelled more accurately than non-endemic species; this effect may be related to the effects of local range size and niche breadth or may be because the environmental gradients inhabited are incompletely sampled in the case of non-endemics (McPherson and Jetz 2007). Overall, we expect species with smaller range sizes, both on local and regional scales, to be modelled more accurately than species with larger ranges. Tests of the effect of range size on model accuracy may be confounded by statistical artefacts. The AUC statistic is a common measure of the accuracy of species distribution models and has been used in many of the studies reviewed here. However, it may be biassed in favour of species with narrow ranges (Lobo et al. 2008).
Only a few studies have considered the effect of migratory behaviour on the accuracy of distribution models. All such studies have focused on birds, most finding that models for migratory species were poorer than models for non-migratory species (Pearce et al. 2001; McPherson and Jetz 2007), probably because the distributions of migratory species are determined by environmental conditions at very specific times of the year (McPherson and Jetz 2007). Conversely, Stockwell and Peterson (2002) found no difference in model accuracy between migratory and non-migratory species, and Mitchell et al. (2001) found that models for migratory bird species were better than models for resident species. No previous study has compared model accuracy between migratory and non-migratory butterfly species. Pöyry et al. (2008) showed that more mobile butterfly species were better modelled than less mobile species, probably because more mobile species can expand their ranges into uninhabited areas more easily and hence occupy a greater proportion of the suitable habitat than less mobile species (but see Pearce et al. 2001). We expect distribution models to be more accurate for resident butterfly species than for migratory species.
There is evidence that both sample size and prevalence (the relative number of presence and absence records) affect the accuracy of distribution models (Manel et al. 1999; Stockwell and Peterson 2002; Brotons et al. 2004; Luoto et al. 2005; Seoane et al. 2005). Therefore, it is important to account for these factors when comparing model accuracy among species (Karl et al. 2002; Huntley et al. 2004; McPherson et al. 2004). Reported effects of prevalence on model accuracy have been mixed, including both positive and negative relationships (Luoto et al. 2005; Brotons et al. 2004), but we expect model accuracy to increase with sample size.
Some authors have demonstrated evolutionary conservatism of ecological niches (Peterson et al. 1999). Furthermore, there may be a substantial heritability in many of the characteristics of species used to explain differences in model accuracy, particularly range size (Jablonski 1987; Hunt et al. 2005; Beck et al. 2006; but see Webb and Gaston 2003). However, to date, only one study has controlled for phylogeny when investigating differences in distribution model accuracy among species (Pöyry et al. 2008). In this case, incorporating phylogeny did not affect the results, but this may not be true for other taxonomic groups, regions or characteristics of species.
In this study, we test the effect of five characteristics of species (local range size, global range size, migratory behaviour, host-plant specialisation and habitat tolerance) on the accuracy of distribution models for butterflies in Egypt, controlling for the potentially confounding effects of sample size and prevalence on model accuracy. Two separate measures of model accuracy were used to minimise the impact of statistical bias on our conclusions. We also control for the influence of species relatedness using independent contrasts.
Materials and methods
Data describing the occurrence of 59 butterfly species in Egypt were compiled as part of Egypt’s BioMAP project (Gilbert and Zalat 2007). The species data were collected between 1829 and 2006, but most records date from the 20th Century (Newbold et al. 2009). We used five environmental variables as predictors: four principal components, which describe altitude and 19 bioclimatic variables from WorldClim (Hijmans et al. 2005), and a categorical land cover variable based on AVHRR satellite data (Hansen et al. 2000). Land cover is classified into 13 categories globally—needleleaf evergreen forest, broadleaf evergreen forest, needleleaf deciduous forest, broadleaf deciduous forest, mixed forest, woodland, wooded grassland, closed shrubland, open shrubland, grassland, cropland, bare ground, and urban areas (Hansen et al. 2000)—eight of which are found in Egypt. All variables had a resolution of 30 arc seconds (~1 km).
Models were built with Maxent Version 2.3 (Phillips et al. 2006). We generated distribution models for 40 species with at least eight occurrence records. The datasets for each species were divided randomly, with half the records used for model building and half for model evaluation. For each species we built ten models, using different subsets of the data to build and evaluate each.
We initially evaluated the models using the area under the receiver operating characteristic curve (AUC). AUC values were calculated using the trapezoid method (Pearce and Ferrier 2000a). Its calculation requires both presence and absence records. We generated 2,500 pseudo-absences, randomly situated in grid cells with no recorded occurrence of a given species. The AUC statistic may be sensitive to the extent of the study area and the proportion of this area that the species inhabits. As an additional evaluation of model performance, we fitted a generalised linear model (GLM) with binomial errors using the presences and pseudo-absences reserved for model evaluation as the dependent variable and the model predicted probability of occurrence at these sites as a single independent variable. The deviance explained by this model was used as a measure of model accuracy. Negative values indicate negative relationships between model predicted probability and observed occurrence. AUC values and deviances explained were averaged across all ten model runs for each species.
We considered six characteristics of species that may affect the accuracy of distribution models: (1) the mean number of presence records used in the ten model runs; (2) whether the species is a migrant, partial migrant or resident in Egypt; (3) whether the species is a specialist or generalist in terms of the host plants it uses; (4) the species’ inhabited range size within Egypt; (5) its global range size (endemic, near-endemic, restricted-range, narrowly distributed or widespread); and (6) its habitat tolerance. Migratory behaviour data were taken from Gilbert and Zalat (2007). Species were defined as specialists if their host plants are confined to one genus and as generalists otherwise, according to Gilbert and Zalat (2007). Maxent outputs a cumulative predicted probability of occurrence for each model between 0 and 100. The mean proportion of grid cells with a predicted value of greater than 50 (averaged across the ten model runs for each species) was used as an index of range size within Egypt. Global range size followed the classifications used in Gilbert and Zalat (2007). The breadth of a species’ habitat tolerance was estimated as the number of land cover categories into which recorded species occurrences fell into.
The results of cross-species comparisons may be confounded by an effect of species relatedness on their niches and on the species characteristics considered. To control for this we calculated independent contrasts for both measures of model accuracy and all six characteristics of species (Harvey and Pagel 1991). One species characteristic (migratory behaviour) was a categorical variable with more than two categories, so we reclassified it into a series of binary variables, one for each category of the original variable. A phylogenetic topology was generated based on published studies (Pieridae: Pollock et al. 1998; Braby et al. 2006; Lycaenidae: Pierce et al. 2002; Pech et al. 2004; Nymphalidae: Brower 2000; Wahlberg et al. 2003; Freitas and Brown 2004; All groups: García-Barros 2000; Wahlberg et al. 2005). In the absence of data describing branch lengths, all branches were assigned a length of one, assuming punctuational evolution (Bro-Jørgensen 2007). We inserted small branches of length 0.0001 in to polytomous clades. The phylogenetic tree was constructed in TreeView 1.6.6 (Page 1996) and modified using Mesquite 1.12 (Maddison and Maddison 2007). The independent contrasts were calculated using Compare Version 4.6b (Martins 2004).
Statistical analysis
We arc-sin transformed model-accuracy measures (both AUC values and deviance explained) to meet assumptions of normality where appropriate. The effects of species characteristics on model accuracy were assessed using analyses of covariance; these analyses took AUC values and deviances explained by the models as the dependent variables, respectively. The six species characteristics were considered as independent variables. Preliminary analysis suggested that two-way interactions did not have a significant effect on model accuracy, so these were excluded from the final analyses.
We used a model selection method based on the approach recommended by Burnham and Anderson (2002). First, we built a global model with all six terms, and candidate models with every combination of terms. AIC scores were extracted for each model and the difference between a model’s AIC value and the lowest value of all models (the AIC difference, Δi) was calculated. Model weight was calculated using the following formula (Burnham and Anderson 2002):
where Δi is the AIC difference of the model in question and Δrs are the AIC differences of the other models. The relative importance of each variable was assessed by summing the AIC weights of all candidate models containing it (Burnham and Anderson 2002), hereafter referred to as the ‘sum of AIC weights’. To test the effect of including species with very small numbers of presence records on the conclusions drawn, we repeated the same analyses considering only the 22 species with more than 20 unique presence records.
Relationships among independent contrasts for model accuracy measures and species characteristics were analysed using Pearson’s correlation tests.
All statistical tests were carried out in SPSS Version 15.0 and R Version 2.6.1 (R Development Core Team 2004).
Results
Models were generally accurate, attaining a mean AUC value of 0.83 ± 0.015 and explained a mean deviance in species occurrence of 23.31 ± 2.98. Predicted range size within Egypt had a strong negative effect on model performance, using both AUC values (sum of AIC weights = 0.921; Table 1; Fig. 1a) and deviances explained by the models (sum of AIC weights = 0.997; Table 2; Fig. 1b) as measures of model accuracy. World range also had a strong negative effect on model accuracy, measured using both AUC values (sum of AIC weights = 0.733; Table 1; Fig. 1c) and the deviance explained by the models (sum of AIC weights = 0.988; Table 2; Fig. 1d). World range and range within Egypt did not correlate significantly with one another (Spearman rank correlation: rs = 0.120, n = 40, P > 0.05). There was limited support for an effect on model accuracy of the number of presence records used to build models (sum of AIC weights = 0.347 and 0.453, for AUC values and deviances explained by models, respectively), migratory behaviour (sum of AIC weights = 0.382 and 0.224), host-plant specificity (sum of AIC weights = 0.316 and 0.318) or habitat tolerance (sum of AIC weights = 0.302 and 0.411). Considering only species with more than 20 unique presence records did not qualitatively alter the results, although migratory behaviour appeared to be a more important determinant of model accuracy in these analyses (Tables 3, 4 in Appendix).
Results of a set of general linear models testing the effect of species characteristics on the accuracy of species distribution models for 40 Egyptian butterfly species, measured using the AUC statistic
Characteristics tested were: the number of presence records used to build models (P), migratory behaviour (M), host-plant specificity (S), predicted range size in Egypt (R), world range size (W) and habitat tolerance (H). Candidate models were built with every possible combination of terms. These models were compared using the approach recommended by Burnham and Anderson (2002), by calculating AIC values for each model, the difference between the AIC for a model and the minimum AIC for all models (Δi), and model weights based on these values. We only present the best models (Δi ≤ 2) here
Effect of species characteristics on model accuracy. Accuracy was measured using both the AUC statistic and the deviance explained by the models. Species characteristics that had a significant effect on model accuracy are shown: predicted range within Egypt (a, b) and world range (c, d). With the exception of the effect of predicted range size within Egypt on AUC scores, all relationships remained significant after accounting for species relatedness using independent contrasts
Results of a set of general linear models testing the effect of species characteristics on the accuracy of species distribution models for 40 Egyptian butterfly species, measured as the deviance explained by the models
Characteristics tested were: the number of presence records used to build models (P), migratory behaviour (M), host-plant specificity (S), predicted range size in Egypt (R), world range size (W) and habitat tolerance (H). Candidate models were built with every possible combination of terms. These models were compared using the approach recommended by Burnham and Anderson (2002), by calculating AIC values for each model, the difference between the AIC for a model and the minimum AIC for all models (Δi), and model weights based on these values. We only present the best models (Δi ≤ 2) here
When species relatedness was accounted for using independent contrasts, world range still showed a significant negative relationship with model accuracy, estimated using both AUC (rp = −0.323, N = 39, P = 0.045) and deviance explained by the models (rp = −0.478, N = 39, P = 0.002). Predicted range within Egypt showed a significant negative relationship with deviance explained by the models (rp = −0.394, N = 39, P = 0.013), but not with average AUC score (rp = −0.110, N = 39, P = 0.506). All other characteristics tested did not have a significant effect on model accuracy after accounting for the relatedness of species (−0.241 ≤ rp ≤ 0.172, N = 39, P > 0.05).
Discussion
Our results confirm that characteristics of species can strongly affect model accuracy, although the factors considered explained a relatively small proportion of the variation in accuracy measures. Of the six characteristics that we tested, two had consistent significant effects on model performance. Disentangling causal mechanisms for patterns such as these is difficult because range size shows relationships with abundance and occupancy (Gaston et al. 2000), and also with ecological characteristics of species, such as dispersal ability and niche breadth (Beck and Kitching 2007). However, our results are consistent with hypothesised relationships between range size and the accuracy of distribution models.
Species with large local range sizes had less accurate models than those with small range sizes. This is consistent with the results of other studies (Stockwell and Peterson 2002; Brotons et al. 2004; Segurado and Araújo 2004; Hernandez et al. 2006). Species with small ranges included both desert species and species inhabiting the Nile Valley and Delta, thus the effect of range size was not an artefact of certain habitats containing better-modelled species. Some authors have suggested that species with smaller ranges have better-defined habitat requirements, making them easier to model (Brotons et al. 2004; Hernandez et al. 2006). However, contrary to the findings of other studies (Boone and Krohn 1999; Pearce et al. 2001; Kadmon et al. 2003; Berg et al. 2004), we found no evidence of an effect of habitat tolerance on the accuracy of species distribution models. A similar study to our own, comparing model accuracy among butterfly species in a temperate environment (Pöyry et al. 2008), also found no effect of niche breadth. Therefore, it would seem that other characteristics of butterfly species may be more important in determining the accuracy of butterfly distribution models, or that aspects of niche breadth that determine model accuracy are not captured by the measures used.
It has been suggested that the AUC statistic may be biased in favour of species that occupy a small proportion of the study area (Lobo et al. 2008), which may explain the existence of negative relationships between range size and model accuracy. However, in our study, the effect of range size was the same for two independent measures of model accuracy, suggesting that the relationship was not an artefact associated with use of the AUC statistic. More generally, the use of pseudo-absences may affect measures of model accuracy (VanDer Wal et al. 2009) and thus relationships between range size and model accuracy. However, we found a strong effect of both global and local range size on model accuracy. While the effect of local range size may be affected by statistical artefacts, the effect of global range size should not.
Species with larger ranges may be modelled less accurately because the study area contains discrete populations that show different responses to the environment (Stockwell and Peterson 2002; Brotons et al. 2004). Although some studies suggest that niches are highly evolutionarily conserved (Peterson et al. 1999), others have found that organisms can adapt their niches very rapidly in certain situations (Knouft et al. 2006). The existence of different populations within species that respond differently to the environment is certainly possible in our study; at least two butterfly species (Carcharodus stauderi and Spialia doris) are known to be represented by two sub-species in Egypt (Gilbert and Zalat 2007). Furthermore, the Nile River, Suez Canal and the mountains of the Eastern and Sinai Deserts may present dispersal barriers for some species, causing isolation of populations.
Global range size also had a strong effect on the accuracy of our models. Predictions for endemic, near-endemic and restricted-range species were better than those for more widespread species. This has been shown before for birds (McPherson and Jetz 2007), and recently for butterflies (Marmion et al. 2008). It has been suggested that endemic species are modelled more accurately because the environmental gradients that they inhabit have been completely sampled, whereas only part of the total inhabited environmental space is sampled for non-endemics (McPherson and Jetz 2007). Alternatively, the effect of global range may be brought about by similar mechanisms to the effect of local range size, i.e. larger-ranged species having locally-adapted populations (e.g. Stockwell and Peterson 2002) or having broader habitat requirements that are more difficult to model (e.g. Hernandez et al. 2006).
Previous studies have suggested that the distributions of specialist species are better modelled than those of generalist species (Hepinstall et al. 2002; Segurado and Araújo 2004; Elith et al. 2006). Ours is the first study to test for this effect in butterflies, and we find little evidence that specialists and generalists differ in the accuracy of their distribution models. Butterflies are dependent on certain plant species as host plants and the distribution of these plants can strongly affect the distribution of the butterflies (Araújo and Luoto 2007, but see Quinn et al. 1998). Therefore, it may be the identity, rather than the number, of host plants that affects the accuracy of butterfly distribution models. If the distribution of a butterfly’s host plant is largely determined by climate and habitat, then we might expect that a model for the butterfly that is based on climate and habitat variables will be more accurate than if the host-plant’s distribution is determined by other factors.
Few studies have considered the effect of migratory behaviour on the accuracy of species distribution models and these have focused on bird species, generally finding that migrant species are modelled less accurately than resident species (Pearce et al. 2001; McPherson and Jetz 2007). If anything, partial migrants had the least accurate models in this study. One possible explanation is that the distributions of both residents and migrants are strongly determined by environmental variables, but that each responds slightly differently to those variables. If partially-migratory species consist of separate populations of migrants and residents, then their distribution models will be less accurate than species that are entirely migratory or entirely resident and respond consistently to the environmental variables. Given the weak trend suggested in our data, more work is needed to explore this phenomenon further.
Several authors have reported a significant effect of sample size on model accuracy (Pearce and Ferrier 2000b; Stockwell and Peterson 2002; Phillips et al. 2004; Hernandez et al. 2006), although this effect has been shown to vary among modelling techniques. In this study we used Maxent to build our models and found no relationship between sample size and model performance. This supports the results of other studies that have shown that Maxent is generally robust to variation in sample size and that it produces accurate predictions even with very small samples (e.g. Hernandez et al. 2008). Most studies of the effects of sample size on model performance (Pearce and Ferrier 2000b; Stockwell and Peterson 2002; Phillips et al. 2004; Hernandez et al. 2006) have experimentally altered sample sizes for one species. We tested the effect of the available sample size across many species. It may be that the completeness of sampling with respect to the environmental gradients rather than sample size alone is most important in determining model accuracy, although Kadmon et al. (2003) found that distribution-model accuracy decreased with the completeness of sampling with respect to climatic gradients.
It is important to account for the effect of species relatedness in comparisons of models across species; otherwise, false conclusions may be drawn regarding the effect of some species characteristics on model accuracy, as is the case in other comparative studies (e.g. Harvey and Pagel 1991). Although accounting for species relatedness had no effect on the conclusions of this study, species distributions, and also some of the species characteristics tested, are known to be evolutionarily conserved (Jablonski 1987; Peterson et al. 1999; Hunt et al. 2005).
The results have important consequences, both for species distribution modelling itself and for conservation biology more generally. It is important to understand why models for different species perform differently before using them to make conservation decisions. This is the first test of differences in accuracy among distribution models of butterflies in an arid environment. The results are generally consistent with those of similar studies of butterflies in other parts of the world, although we present the first test of the effects of migratory behaviour and host-plant specialism on the accuracy of models for butterfly species. It is important to note that the factors that determine species distributions vary according to the scale of analysis (Whittaker et al. 2001), and thus the characteristics of species that affect distribution-model accuracy may also differ. We also emphasise that it is important to control for phylogeny when conducting cross-species comparisons like this one. Although there was substantial variation among species in model accuracy, accurate models were produced for many species, confirming the value of such models in conservation ecology.
Acknowledgments
We thank Italian Cooperation (Debt Swap) for funding the BioMAP Project, Dr. Mustafa Fouda (Director of the Nature Conservation Sector, EEAA) for facilities and comments on the work, all the BioMAP staff (Ahmed Yakoub, Alaa Awad, Muhammed Sherif, Shama Omran, Shaimaa Esa, Yasmin Safwat, Nahla Ahmed, Esraa Sabre), Dr. Abd El Aal Attia for help during dataset preparation and preliminary analysis. Two anonymous reviewers made valuable comments on an earlier draft of this paper. This work was supported by the Natural Environment Research Council (grant number NER/S/A/2006/14170).
To check that our conclusions were not biased by including species with very small numbers of presence records, we repeated the analyses of the effect of characteristics of species on distribution model accuracy, considering only species with at least 20 unique presence records. See Tables 3 and 4.
Results of a set of general linear models testing the effect of species characteristics on the accuracy of species distribution models for 22 Egyptian butterfly species with at least 20 unique presence records, measured using the AUC statistic
Characteristics tested were: the number of presence records used to build models (P), migratory behaviour (M), host-plant specificity (S), predicted range size in Egypt (R), world range size (W) and habitat tolerance (H). Candidate models were built with every possible combination of terms. These models were compared using the approach recommended by Burnham and Anderson (2002), by calculating AIC values for each model, the difference between the AIC for a model and the minimum AIC for all models (Δi), and model weights based on these values. We only present the best models (Δi ≤ 2) here
Results of a set of general linear models testing the effect of species characteristics on the accuracy of species distribution models for 22 Egyptian butterfly species with at least 20 unique presence records, measured as the deviance explained by the models
Characteristics tested were: the number of presence records used to build models (P), migratory behaviour (M), host-plant specificity (S), predicted range size in Egypt (R), world range size (W) and habitat tolerance (H). Candidate models were built with every possible combination of terms. These models were compared using the approach recommended by Burnham and Anderson (2002), by calculating AIC values for each model, the difference between the AIC for a model and the minimum AIC for all models (Δi), and model weights based on these values. We only present the best models (Δi ≤ 2) here

Tim Newbold, Francis Gilbert, Samy Zalat, Ahmed El-Gabbas, Tom Reader 2009
Climate-based models of spatial patterns of species richness in Egypt’s butterfly and mammal fauna
Journal of Biogeography 36: 2085-95
10.1111/j.1365-2699.2009.02140.x
Abstract
Aim: Identifying areas of high species richness is an important goal of conservation biogeography. In this study we compared alternative methods for generating climate-based estimates of spatial patterns of butterfly and mammal species richness.
Location: Egypt.
Methods: Data on the occurrence of butterflies and mammals in Egypt were taken from an electronic database compiled from museum records and the literature. Using Maxent, species distribution models were built with these data and with variables describing climate and habitat. Species richness predictions were made by summing distribution models for individual species and by modelling observed species richness directly using the same environmental variables.
Results: Estimates of species richness from both methods correlated positively with each other and with observed species richness. Protected areas had higher species richness (both predicted and actual) than unprotected areas.
Main conclusions: Our results suggest that climate-based models of species richness could provide a rapid method for selecting potential areas for protection and thus have important implications for biodiversity conservation.
Introduction
In 2002, the Convention on Biological Diversity adopted the target of significantly reducing the rate of biodiversity loss by 2010. Monitoring biodiversity trends requires knowledge of the distributions of species and of patterns in species richness. However, such data are limited, especially in the tropics (Anderson et al., 2003) and in arid regions (Stockwell & Peters, 1999).
Using statistical methods to predict the distribution of species is an approach that shows great promise (Wintle et al., 2005). Distribution models relate known species occurrence data to environmental variables, such as climatic variables and habitat types. Both traditional statistical techniques and newer methods have been shown to model individual species distributions with a high degree of accuracy (Elith et al., 2006). Some authors have experimented with summing individual predictions of species distributions to estimate species richness. For example, Garcia (2006) produced distribution models for 267 reptile and amphibian species in Mexico and summed them to produce a prediction of species richness. However, when large numbers of species are involved, summing distribution models to generate species richness predictions may be time-consuming (Gioia & Pigott, 2000). An alternative approach is to model species richness directly. Whilst there have been many attempts to find climatic and habitat-related determinants of species richness patterns (e.g. Kivinen et al., 2007; Levinsky et al., 2007), no studies have explicitly compared summed distribution predictions with models of species richness per se, although Gioia & Pigott (2000) used both approaches. Such a comparison will be very useful for conservation biologists attempting to understand spatial patterns of biodiversity.
Many studies have investigated correlates of species richness, often finding climate variables to be good correlates of observed patterns (e.g. Hawkins et al., 2003). Arid environments are under-studied in this respect (but see van Rensburg et al., 2002; Schmidt et al., 2008). The mechanistic explanation for these relationships remains a matter of debate, and the conclusions of any study of patterns of species richness are strongly affected by the spatial scale at which they are conducted (Field et al., 2009). Currie et al. (2004) explored three hypotheses for climate-based patterns in species richness, concerning energy, the climatic tolerance of species and speciation rates. They did not find unequivocal support for any of these hypotheses in the literature. At broad scales, historical factors (Qian & Ricklefs, 2000) and the distribution of resources (Araújo & Luoto, 2007) may play important roles in determining species richness. At finer scales, competition (Anderson et al., 2002), metapopulation dynamics (Hanski, 1991) and human disturbance (Uehara-Prado et al., 2007) have also been shown to exert a significant influence on species richness.
Several studies have shown that butterfly and mammal species richness correlate with climate and habitat variables in temperate and tropical regions, at both local scales (Turner et al., 1987; Kivinen et al., 2007; Kuussaari et al., 2007) and at broader scales (Nogués-Bravo & Araújo, 2006; Algar et al., 2007; Levinsky et al., 2007). However, to date very few studies have investigated correlates of mammal and butterfly species richness in an arid environment (but for mammals see Andrews & O’Brien, 2000).
In this paper we report the results of a study of butterfly and mammal species richness in Egypt, a country with a typical arid-environment flora and fauna. Egypt has two endemic and two near-endemic butterfly species, and also three endemic subspecies (Larsen, 1990). The mammal fauna includes four endemic and 10 near-endemic species. We sought to identify environmental correlates of species richness at a local scale in an arid environment. We also asked whether estimates of the species richness of Egyptian butterflies and mammals derived from models of species richness had a good match with estimates made by summing individual models of the distribution of species, and whether both these estimates matched observed patterns of species richness.
One application of models of species richness is in assessing the effectiveness of protected areas. Global estimates of the effectiveness of protected areas generally suggest poor coverage of biodiversity (Chape et al., 2005). Country-level studies have often found species richness to be no higher in protected areas than in unprotected areas (e.g. Pawar et al., 2007; Traba et al., 2007; but see, e.g., Lee et al., 2007). Egypt has 27 current or proposed protected areas, covering a total of 11% of its land surface (Egyptian Environmental Affairs Agency, 2007). All these have been gazetted since 1983, mostly at the recommendation of scientists familiar with Egypt’s biodiversity. As such, we may expect them to show better coverage of biodiversity than protected areas in other countries. We tested whether protected areas in Egypt have higher species richness than unprotected areas.
Materials and methods
Species and climate data
Species occurrence data were compiled as part of Egypt’s BioMAP project (see http://www.biomapegypt.org/ for more details). The butterfly dataset consisted of 1898 records for 59 species, mostly from museum specimens and the sparse literature on Egyptian butterflies (Larsen, 1990; Gilbert & Zalat, 2007). Records were made between the years 1829 and 2006, but most were from the 20th century (see Fig. S1 in Appendix S1 in Supporting Information). All extant specimens in Egyptian collections were re-identified according to the latest taxonomic opinion; Larsen (1990) had already reviewed and checked most other records. Coverage of Egypt is patchy, but probably fairly representative in the sense that all the main areas where butterflies occur have been sampled, whilst the huge areas of the Great Sand Sea where no butterflies are thought to occur have not been sampled. The main lacuna in collecting effort is the Qattara Depression in the northern part of the Western Desert.
The mammal data consisted of 4718 records for 103 species, also taken from museum and personal records, unpublished reports and the published literature. Mammal records were made between the years 1580 and 2007, but the majority fell in the second half of the 20th century (see Fig. S1 in Appendix S1). The identification of every record and Egyptian specimen was checked or rechecked, according to the latest taxonomic opinion (Wilson & Reeder, 2005), by the top mammal taxonomist in Egypt (Dr M. Bassiouny, Al Azhar University, Cairo). Coverage of Egypt is very good, due to systematic and extensive collecting in the period 1950–80 coupled with careful evaluation of taxonomy (see Osborn & Helmy, 1980).
All locations were mapped as accurately as possible using a bespoke gazetteer developed by the BioMAP project over 3 years of collating Egyptian biodiversity records. The maximum error of each location was estimated by the point radius method (see Wieczorek et al., 2004) and excessively inaccurate records rejected.
Climatic predictors were taken from the WorldClim version 1.4 dataset at a resolution of 30 arcsec (c. 1 km) (Hijmans et al., 2005). This source contains 20 variables describing aspects of temperature, precipitation and elevation. We also used a new Egyptian geological habitat map (hereafter referred to simply as habitat) (A. A. Hassan, unpublished data). Habitat was classified into 11 classes (sea, littoral coastal land, cultivated land, sand dune, wadi, metamorphic rock, igneous rock, gravels, serir sand sheets, sabkhas and sedimentary rocks) based on remote sensing and extensive ground truthing. In preliminary analyses we also experimented with topographical predictors (slope and aspect). However, these variables did not significantly improve model accuracy and were excluded from the final analyses.
To ensure that the species records were not environmentally biased, we tested their coverage of the main environmental gradients. Butterfly and mammal sites were treated separately for this analysis. Four principal component variables based on the 20 climatic and elevation variables were each classified into four categories using Jenks (1967) natural breaks, giving 256 unique combinations of categories (henceforth called areas of climatic space). The number of areas of climatic space that contained at least one survey location was calculated. To estimate the number of areas that would have been visited if sampling was completely random, we generated 100 sets of random survey points of the same number as the real sampling locations and calculated the number of areas represented. To evaluate the coverage of habitat classes by sampled locations, a chi-square test was performed to assess whether sites fell into significantly (α = 0.05) different proportions of classes than expected by chance.
Modelling species richness
We modelled the species richness of butterflies and mammals separately, using two methods. First, we summed predictions of the distribution of individual species, using a resolution of 30 arcsec. We made initial distribution models for the 40 butterfly species and 68 mammal species with at least eight records of occurrence, using Maxent version 2.3 (Phillips et al., 2006). Maxent is designed for use with species datasets that contain only records of presence (Phillips et al., 2006) and thus may be particularly useful when species data are taken from museum collections. Maxent attempts to find a probability distribution that is as close to uniform as possible with the constraint that the expected value of each environmental variable (the sum, across all grid cells, of the product of the probability of occurrence and the value of the environmental variable) must equal the mean value at the presence points (the empirical average). Regularization relaxes this constraint such that the expected value of each environmental variable may fall within a predefined margin around the empirical average, preventing over-fitting of the model (Dudík et al., 2004). The algorithm runs until improvement in model accuracy at each iteration falls below a set threshold (convergence) or until a maximum number of iterations have been performed. Full details are given in Phillips et al. (2006). We used the 19 climatic variables, elevation and habitat as predictor variables. Linear and quadratic terms were fitted for continuous variables. We used default values for all parameters (a regularization value of 1, a convergence threshold of 0.00001, a maximum of 500 iterations and a sample of 10,000 points to characterize the background environment). Ten initial models were made for each species. For each model, the species data were randomly divided into half for building the model and half for evaluating the model. The accuracy of each model was assessed using the ‘area under the receiver operating characteristic curve’ (AUC) statistic, as calculated within the Maxent procedure. Following the recommendations made by Swets (1986), we eliminated five butterfly species and seven mammal species with mean AUC scores of < 0.7. This left 35 butterfly species, including one of the two endemic species and both near-endemic species, and 61 mammal species, including three of the four endemic species and five out of 10 near-endemic species.
A single final model was then made for each of the remaining species, again at a resolution of 30 arcsec, using the same modelling protocol, but this time using all presence records. The output of statistical models varies among species according to the relative numbers of presences and absences in the species data (prevalence) (Manel et al., 2001). Therefore, simply summing the output of individual distribution models may bias estimates of species richness in favour of taxa with many records. It is better to convert the model output into a binary prediction of presence or absence around a threshold value. Many methods have been proposed for choosing appropriate thresholds. For datasets consisting only of presences, Pearson et al. (2004) recommended using a threshold that maximizes sensitivity (the percentage of presences correctly predicted as being present at a given threshold). Here we used a threshold that resulted in a sensitivity of 95%. Once the models had been converted to binary predictions of presence or absence, they were summed across all species to give an estimate of species richness.
The second method of modelling species richness was to model observed species richness values directly. This part of the study was concerned with the total number of species recorded in each cell rather than individual records of species. Therefore, we used a resolution of 0.5°, because at the finer resolution most cells contained no records. Observed species richness was calculated from the original survey data in diva-gis 5.2 (http://www.diva-gis.org/). A species was considered present in a cell if it had been recorded at least once. Species richness was modelled using generalized linear models (GLMs) with the same independent variables as in the species distribution models. Following an inspection of the residuals of a general linear model and consideration of dispersion, the most appropriate family of GLMs was chosen for each model. The variables were resampled to the coarser resolution using bilinear interpolation. In bilinear interpolation, the values of the four nearest grid cells to the target cell are averaged after being weighted according to their distance to the target cell. Fitting too many independent variables in GLMs may result in overfitting and the selection of nonsense variables in the final model (Wintle et al., 2005). To avoid these problems, we performed a principal components analysis (PCA) on the 19 climatic variables and elevation across all 406 of the 0.5° cells. Components with an eigenvalue > 1.0 were retained as new predictor variables. In the GLMs, linear and quadratic terms were fitted for each of these components. We constructed two separate models of species richness: (1) assuming that cells with no species records had a species richness of zero; and (2) excluding cells with a recorded species richness of zero. In the first we fitted recorded species richness values of all 406 of the 0.5° grid cells in the study area. We used a GLM with negative binomial errors (NB-GLM) and the log link (Crawley, 2002; Venables & Ripley, 2002). For the second model, since some cells with a recorded richness of zero may occur simply because they have not been sampled and the results could be biased by the inclusion of false zero values, we fitted the species richness values of 0.5° cells with at least one record of the taxonomic group in question – 100 cells for butterflies and 196 cells for mammals. A GLM with Poisson errors (P-GLM) and the log link (McCullagh & Nelder, 1989) was used.
Statistical analysis
The agreement between fitted values of species richness generated using the different methods was tested using Spearman’s rank correlation tests. For comparison, the species richness prediction generated by summing individual distribution models was resampled from its original resolution of 30 arcsec to a resolution of 0.5° using bilinear interpolation. Thus, all tests compared species richness across all 362 of the 0.5° grid cells with an estimate of species richness by all three models. These comparisons included cells with no species records; these cells were assumed to have a species richness of zero. We also repeated the same correlation tests using only cells that had at least one record of a species in the group being considered.
We tested whether Egypt’s protected areas network represented butterfly and mammal species richness well by comparing estimated (using the distribution model-sum method) and observed species richness inside and outside protected areas at 2000 points, randomly situated in 1-km cells throughout the study area. These points were generated using Hawth’s analysis tools for ArcMap 9.1 (Beyer, 2004). We also compared both estimated and observed richness of endemic and near-endemic species inside and outside protected areas. For this comparison, we grouped mammals and butterflies because the number of endemic species was small.
The NB-GLMs and P-GLMs were built using the glm (Poisson errors) and glm.nb (negative binomial errors) packages in R 2.6.1 (R Development Core Team, 2004). For both, a manual backward stepwise selection procedure was used to remove terms that did not significantly improve the deviance explained, until a minimum adequate model was obtained (Crawley, 2002). All other analyses were carried out using spss 15.0 (SPSS Inc., Chicago, USA). The comparison of actual and predicted species richness inside and outside protected areas was undertaken using a Mann–Whitney U-test.
Results
The species occurrence data showed no biases in environmental space that would preclude accurate modelling. Sampled sites covered the full range of values of each of the principal components based on the climatic variables. The butterfly sampling locations fell into 44 of 256 areas of climatic space, 84.1% of the number expected by chance. Sites that were sampled for mammals covered 76 areas of climatic space, 107.5% of the number expected by chance. Sites for both butterflies and mammals fell into significantly different proportions of habitat types than expected by chance (χ2 = 1035, d.f. = 9, p 0.001 and χ2 = 2248, d.f. = 9, p 0.001, respectively). For both butterflies and mammals, littoral coastal areas, cultivated land, wadis (dry desert valleys), areas of metamorphic rock and areas of igneous rock were sampled more often than expected by chance. Sand dunes, gravels, serir sand sheets and areas of sedimentary rock were sampled less often than expected by chance. Sixty-three of 333 butterfly sites and 200 of 1395 mammal sites fell inside Egypt’s protected areas. A map of Egypt’s protected areas and the sites that were sampled for mammals and butterflies is given in Fig. 1(a).
(a) Sites where mammals (circles) and butterflies (triangles) were sampled, and the location of Egypt’s protected areas (grey shading). Maps of (b) predicted butterfly and (c) mammal species richness generated by summing individual predictions of the distributions of species. Lighter tones indicate high predicted species richness and darker tones indicate lower species richness. The distribution predictions were made using Maxent.
The final distribution models (those using all the species occurrence data) achieved mean AUC values between 0.863 and 0.999 for butterfly species and between 0.831 and 0.999 for mammal species. The average relative contribution of habitat, elevation and the 19 climatic variables to the final distribution models is shown in Fig. S2 in Appendix S1, and the contributions of variables in the models for each species are given in Appendix S2. Habitat and elevation were important in explaining the distributions of both butterflies and mammals. Among the climatic predictors, temperature-related variables were more important than precipitation-related variables in determining butterfly distributions, while for mammals, annual and maximum precipitation variables were also important.
The predictions of species richness made using the first method (summing distribution models for individual species) are mapped in Fig. 1(b,c). The models of species richness generated using this method correlated positively and significantly with observed species richness (Table 1, Fig. 2a). The second method of estimating spatial patterns of species richness was to model species richness values directly. The PCA of the 19 climatic variables and elevation produced four components with eigenvalues > 1.0. All original climate variables were represented in at least one of the extracted components (Table S1 in Appendix S1). Scores on the first principal component increased with increasing maximum temperature and decreasing precipitation annually and at otherwise wet times of year (Table S1 in Appendix S1). The second component increased with increasing annual temperature and increasing temperature during cooler periods of the year. The third component described increasing elevation, decreasing annual temperature and increasing precipitation at drier times of the year. The fourth component increased with decreasing temperature during dry periods, increasing precipitation (annually and during cold times of the year) and decreasing minimum precipitation. The models fitting species richness values for all 406 of the 0.5° cells (NB-GLM), which included cells with a recorded species richness of zero, explained 16.3% of the deviance in the species richness of butterfly species and 21.3% of the deviance in mammal species richness. For butterflies, only the linear term of PC1 (describing mainly precipitation but also maximum temperature) and habitat had a significant effect on species richness (Table 2). For mammals, habitat, the quadratic term of PC1 and both terms of PC2 (describing several temperature-related variables) and PC4 (describing variables related to extremes of temperature and rainfall) had a significant effect on species richness (Table 2). Estimated species richness according to the NB-GLMs correlated significantly and positively with observed species richness (Table 1, Fig. 2b). The models fitting species richness values only for 0.5° cells with at least one species record (P-GLM) explained 19.1% of the deviance in butterfly species richness and 18.3% of the deviance in mammal richness. For butterflies, both terms of PC1 and PC4, the quadratic term of PC3 (high values of which indicate high-altitude areas with rainfall all year round) and habitat were all significantly related to species richness (Table 3). For mammals, both terms of PC1, the quadratic terms of PC2 and PC4, the linear term of PC3 and habitat were significant predictors of species richness (Table 3). Species richness estimates from these models also correlated significantly and positively with observed species richness (Table 1, Fig. 2c). Across all 0.5° grid cells, the estimates made using the different modelling methods correlated significantly with each other (Table 1, Fig. 3).
Table 1. Correlations among fitted values of each of the three models of species richness of butterflies (B) and mammals (M) in Egypt, and correlations between these fitted values and observed species richness. The three models of species richness were: (1) summed distributions – distribution models were built for each species at 30 arcsec resolution using Maxent, then summed to estimate species richness; (2) NB-GLM – species richness values of all 0.5° cells were fitted using a generalized linear model (GLM) with negative binomial errors; (3) P-GLM – species richness values of sampled cells were fitted using a GLM with Poisson errors. Correlations were calculated both for all cells and for sampled cells only. Species richness values cannot be considered independent in the presence of spatial autocorrelation; the effective sample size is reduced in the presence of such non-independence. Therefore, the minimum sample sizes at which the reported correlation coefficients would remain significant (at α = 0.05) are given in brackets after the correlation coefficient.
Correlations between observed species richness and predicted species richness of butterflies and mammals in Egypt estimated using each of the three models: (a) predictions of the distribution of each species, produced using Maxent, were summed; (b) recorded species richness values of all grid cells were modelled using a generalized linear model with negative binomial errors; and (c) recorded species richness values of sampled grid cells were modelled using a generalized linear model with Poisson errors. Table 2. Results of generalized linear models (GLMs) with negative binomial errors, fitting the observed species richness of butterflies and mammals in Egypt of all 0.5° grid squares as the dependent variable, with habitat and four bioclimatic principal component axes (linear and quadratic terms) as independent variables. Significant terms are shown in bold. LR, likelihood ratio.
Table 3. Results of generalized linear models (GLMs) with Poisson errors, fitting the observed species richness of butterflies and mammals in Egypt of sampled 0.5° grid cells only (i.e. excluding zero values), with habitat and four bioclimatic principal component axes (linear and quadratic terms) as independent variables. Significant terms are shown in bold.
Correlations between the different models used to predict species richness patterns of butterflies and mammals in Egypt: (a) between the sum of individual species distribution models and the generalized linear model with negative binomial errors (NB-GLM) of species richness values; (b) between the sum of individual species distribution models and the generalized linear model with Poisson errors (P-GLM) of species richness values; and (c) between the NB-GLM and P-GLM models of species richness values. y = x lines are shown for reference.
Observed mammal and butterfly species richness values correlated significantly with each other at sites where at least one butterfly species and one mammal species had been recorded (rs = 0.615, n = 97, p 0.001). Predicted species richness (estimated using the distribution model-sum method) also correlated strongly and significantly between butterflies and mammals (rs = 0.920, n = 362, p 0.001).
Across a random sample of 2000 1-km grid cells, predicted species richness, estimated by summing individual modelled species distributions, of both butterflies (Mann–Whitney test: U = 76,100, n = 1995, p 0.001) and mammals (Mann–Whitney test: U = 70,300, n = 1995, p 0.001) was significantly higher inside protected areas than outside (Fig. 4a). The observed species richness was also significantly higher inside protected areas than outside for both butterflies (Mann–Whitney test: U = 111,000, n = 1995, P = 0.016) and mammals (Mann–Whitney test: U = 80,700, n = 1995, p 0.001) (Fig. 4b). Predicted (Mann–Whitney test: U = 105,000, n = 1963, P = 0.028) and observed (Mann–Whitney test: U = 102,000, n = 1963, P = 0.001) richness of endemic and near-endemic species (mammals and butterflies combined) was significantly higher inside protected areas than outside.
(a) Comparison of predicted species richness (mean ± SE) of butterflies and mammals in Egypt, estimated by summing individual species distribution models, between protected areas and unprotected areas. (b) Comparison of observed species richness (mean ± SE) between protected areas and unprotected areas.
Discussion
We found significant relationships between species occurrence, species richness and the climate and habitat variables that we used. Beale et al. (2008) found that relationships between species occurrence and environmental variables were no better than expected by chance. However, we found that habitat and climate variables both had a significant effect on butterfly and mammal distributions, and on patterns of species richness, at least at the local scale at which we carried out this study. This finding is consistent with previous studies of butterflies and mammals, where climate and habitat have been identified as good predictors of richness, both at continental and at local scales (e.g. Nogués-Bravo & Araújo, 2006; Algar et al., 2007; Kivinen et al., 2007; Levinsky et al., 2007). The association with habitat may reflect the effect of variation in plant communities on animal species distributions. Butterflies and herbivorous mammals are directly dependent on plants for food, whilst other mammal species may rely on certain vegetation types indirectly, for example through the availability of herbivorous prey. Temperature variables appear to have a particularly strong effect on butterfly species, although causality cannot be inferred from correlative models. Similar relationships have been noted before (Turner et al., 1987) and could be brought about by direct effects of temperature on thermoregulation, or indirectly through climate-driven variation in habitat diversity or plant productivity.
Many other factors in addition to climate may affect species richness, including competition (Anderson et al., 2002), the availability of host plants (Araújo & Luoto, 2007), metapopulation dynamics (Hanski, 1991), human disturbance (Uehara-Prado et al., 2007) and other environmental variables such as soil type (Kuussaari et al., 2007), although some of these factors are likely to play a role in determining species richness only at larger spatial scales than were studied here (Whittaker et al., 2001). Given all these non-climatic determinants of species richness patterns, it is not surprising that only a relatively low proportion of the variation in species richness was explained by the models, and that the correlations between modelled and observed species richness were only moderately strong. Some progress is being made towards including factors other than climatic ones in species distribution models (e.g. Araújo & Luoto, 2007) and this must remain a priority for improving the accuracy of predictions. However, the need to identify areas to conserve is urgent and we cannot wait to act until the most accurate models possible have been built for every species. Climate-based models match observed distributions well and are quick and easy to build for a large number of species.
Another reason for the relatively low explanatory power of the models may be that species inventories in sampled cells were incomplete. This seems likely, given that it may be necessary to visit a site many times before absence can be inferred with confidence (MacKenzie et al., 2002). In the case of the NB-GLMs, the inclusion of cells with no records of species presence may have introduced false absences to the models. This is especially likely for the butterfly models, because surveying was less extensive. Ground-truthing will be required to assess the extent to which mismatches between modelled and observed species richness are due to incomplete inventories or to errors in the models.
Across all grid cells in the study area, the three methods produced models that showed positive correlations with observed species richness and with each other, suggesting that they could all be used to predict the species richness of unknown areas from limited data on the distributions of species, an application that would be of great value for conservation. There were some differences among the models though. Species occurrence and richness data often contain many absences or zero values, especially datasets for small or cryptic species with a low probability of detection (MacKenzie et al., 2002). This can bias the parameter estimates of statistical models (Martin et al., 2005). The weaker correlations between observed species richness and species richness estimated using the NB-GLMs (those that included cells with recorded species richness of zero), especially for butterflies, may be caused by the inclusion of false absences. This conclusion is further supported by the observation that the NB-GLMs produced much lower estimates of species richness than the other two methods (see Fig. 3). The mammal data covered a much larger proportion of both geographical and environmental space than the butterfly data, suggesting that recorded species richness values of zero were more reliable.
Summing the individual distribution models produced the best estimates of species richness. The relationship between predicted and observed species richness across sampled cells was slightly weaker for butterflies, which may be a result of model bias caused by false absence records. This would be concerning, given that Maxent is designed to be used with datasets containing only presences (Phillips et al., 2006) and is a possibility that deserves further attention. However, for both butterflies and mammals, summing distribution models produced the estimates of species richness that matched observed species richness most closely. The P-GLM models, which did not include cells with a recorded species richness of zero, generated more accurate estimates of species richness than the NB-GLM models. GLMs are quick and easy models to build, making them a good first choice for modelling species richness patterns.
Some previous work has indicated good spatial agreement among different groups in their species richness at regional scales (Qian, 2007), although a global study of bird, mammal and amphibian species richness found the converse (e.g. Grenyer et al., 2006). The results of this study show that, at least at a local scale within a single country, butterfly and mammal diversity correlate strongly and positively.
In contrast to the findings for many other countries and taxonomic groups (e.g. Evans et al., 2006; Pawar et al., 2007; Traba et al., 2007), Egypt’s protected areas network seems to be effective in representing butterfly and mammal diversity. In many parts of the world, protected areas have historically included land that has relatively little commercial value; such areas do not necessarily represent the best choice in terms of conserving biodiversity (Margules & Pressey, 2000). Egypt’s protected areas network is relatively new and the areas were chosen with the aid of knowledge about the country’s biodiversity. Given such knowledge, and the ability to overcome conflicting interests over land use, it seems that good coverage of biodiversity can be achieved. On the other hand, large areas of the Nile Valley and Delta were predicted to have relatively high butterfly diversity but are not yet protected, suggesting that although great progress has been made towards protecting species-rich areas, more could still be done. The Nile Valley contains land of high economic value, and setting aside areas to be protected may present a challenge. It is important to note that species richness is only one measure of the importance of conserving different areas. Some authors have suggested using taxonomic uniqueness (e.g. Kershaw et al., 1995), complementarity (Margules & Pressey, 2000) or threat (Wilson et al., 2007) instead. Many of Egypt’s endemic and near-endemic mammal and butterfly species were included in the models and the richness of endemic and near-endemic species was higher inside protected areas than outside, but a more comprehensive assessment of the protected areas should consider a number of different criteria.
In summary, we have shown that seemingly accurate estimates of species richness can be made using patchy, incomplete data, allowing us to predict the species richness of sites that have not been surveyed. The three models of species richness were largely similar, although the model based on individual distribution models produced the most consistently accurate results. A similar comparison of the same three models in different regions and for different species would be useful in establishing the general reliability of the approach. Models based on species richness itself, rather than individual species distributions, may be useful when species identity is unknown, for example when using species richness estimators. The results are important for conservation, given the urgency with which we must identify areas that need to be protected, although similar comparisons of species richness models for more taxonomic groups and for a broader geographical region would be useful.
Acknowledgements
We thank the Italian Cooperation (Debt for Development Swap) for funding the BioMAP Project; Mustafa Fouda (Director of the Nature Conservation Sector, Egyptian Environmental Affairs Agency) for facilities and comments on the work; all the BioMAP staff (Ahmed Yakoub, Alaa Awad, Muhammed Sherif, Shama Omran, Shaimaa Esa, Yasmin Safwat, Nahla Ahmed, Esraa Saber); Mohamed Basuony and Abd El Aal Attia for help during dataset preparation and preliminary analysis; and Richard Field and four anonymous referees for valuable comments on an earlier draft. This work was supported by the Natural Environment Research Council (grant number NER/S/A/2006/14170).
Biosketch
The behavioural ecology group at Nottingham University (see http://www.nottingham.ac.uk/~plztr/groupsite/) studies a wide range of topics including conservation ecology and biodiversity. BioMAP (http://www.biomapegypt.org/) spent 3 years gathering species records from museums, natural history collections and the literature, compiling them into an electronic database. It also promoted wildlife conservation in Egypt. The second part of the BioMAP project will commence this year.

Patrick Brading, Ahmed El-Gabbas, Samy Zalat & Francis Gilbert 2009
Biodiversity economics: the value of pollination services to Egypt
Egyptian Journal of Biology, 11, 45-51
Abstrac: Pollinator populations are under severe pressure worldwide because of man-made intensification in land use, including the use of pesticides and fertilizers. The majority of wild and crop plants are fully or partially dependent on pollinators for their reproduction. Loss of pollinators has already caused measurable declines in the populations of many wild plants in Europe. Many Egyptian crops are fully or partially dependent on pollinators for their yields, and data exist on the market values of Egyptian crops. We therefore use these to estimate the 2004 costs to the Egyptian economy of a catastrophic loss of pollinators. The annual cost to the Egyptian economy of losing its pollinators would be approximately LE 13.5 billion ($2.4 billion), 3.3% of the 2003 GDP.
Introduction: Our lives are heavily dependent upon the planet's biodiversity and the ecological systems that it supports. The many products (e.g. raw materials such as timber) and services (e.g. climate regulation) provided by these ecosystems are not only essential to our own survival but also to the functioning of the Earth's life-support system (Millenium Ecosystem Assessment 2005).
Due to the difficulty of placing a realistic monetary worth on ecosystems, their services are not given adequate importance when making policy decisions (van Jaarsveld et al. 2005). Agriculture in particular has many obvious dependencies on natural services provided by the ecosystem. Ironically, however, agriculture is one of the main driving forces behind the decline of biodiversity (UNEP 2007).
Pollination is a prime example of a supporting service that is being negatively affectedby agricultural practices, as well as by other factors such as global warming and urbanisation (Klein et al. 2007). Pollination is essential to most plants for reproduction, including commercial crops. This ecosystem function is provided by many wild pollinator species. There have been worldwide declines in pollinator diversity (Dias et al. 1999, Klein et al. 2007), with declines identified in at least one region or country on every continent (except Antarctica), including the UK and Netherlands in Europe (Biesmeijer et al. 2006).
The types of pollinators in decline include wild bees (social and solitary), domesticated honeybees, hoverflies, butterflies, bats, hummingbirds and other small mammals. The causes of these declines in pollinator biodiversity are almost certainly related to changes in landuse (Klein et al. 2007).
Pollinators require local floral diversity and nesting sites in order to persist in the unnatural environment of farmed land, but loss of natural habitat (usually related to land use practices) prevents this. Agricultural intensification leads to loss and fragmentation of natural pollinator habitat, while climate change, introduction of alien plants and competition with non-native fauna adds to the pressure placed on pollinator populations.
Without the service provided by pollinators, many plant species would be driven to extinction, and cultivation of many modern crops would be impossible. Many crops are wholly dependent on cross-pollination (such as melons and squash) by pollinators, while other crops show significant yield increases when cross-pollinated instead of self-pollinated (such as Brading et al.: apples, tomatoes and cotton). It has been estimated that pollination is responsible for as much as 30% of agricultural food production (UNEP 2007), and in some cases pollination services may contribute as much or more to yields than fertilisers. Due to its ability to dramatically improve yields, the economic value of natural pollination worldwide is thought to be between US$65 and US$70 billion each year (Dias et al. 1999). Inadequate pollination can not only reduce yields, but may also delay them and be the reason for inferior fruit production.
Domesticated honeybees remain the world’s most important pollinators (Klein et al. 2007), but even they are declining and disappearing for no obvious reason (BBC 2007). Without wild pollinator species, current levels of agricultural productivity are under threat.
The International Pollinators Initiative (Dias et al. 1999) was adopted by the Convention on Biological Diversity because of the perceived threat to such a valuable ecosystem service. However, action by politicians and decision-makers is hampered by a lack of estimates of the true value of this ecosystem service. Several methods have been proposed that try to give a monetary value to ecosystem services, none of them perfect (de Groot et al. 2002, Chee 2004). Direct market valuation is the exchange value that services have in trade.
While straightforward, the value is only what the product is worth to a buyer, and omits other less direct values of the services (de Groot et al. 2002). However, it is simple to understand and clear to apply, especially for pollination where the value of the product is often available.
Stimulated by a recent review on experimental evaluations of the impact of pollination on crops (Klein et al. 2007), we use here the Direct Market Valuation approach to estimate the economic value of pollination to the Egyptian economy.
Materials & Methods
Egypt’s main arable output covers 70 different plants, including non-consumed field crops (such a cotton), fruits, nuts and vegetables. These plants differ in their reliance on pollinators for successful fruit and seed setting, from full dependence (e.g. watermelon, melon, custard apples) to total independence (e.g. date, grape, maize, olive). The review by Klein et al. (2007) places the available information for each crop into one of four categories of the impact of pollinator loss on yield: essential (reduction of >90%), high (40-90% loss), modest (10-40%), little (0-10%) and none (0%). For calculation, we used the midpoint of these ranges: 95, 65, 25, 5 and 0 respectively. Although these are approximations, when summed over all the types of Egyptian crops, the final figure is likely to be a reasonable estimate.
Klein et al. (2007) was based on a worldwide review rather than an Egyptian-specific one, and it would be very useful to have an equivalent review of Egyptian pollination studies:
alas, such a review does not exist. The estimated reductions represent the average loss in yield obtained in all the various experiments carried out on any one crop type anywhere in the world.
The figures might well be different under Egyptian conditions, but until the relevant experiments and review have been done, we do not know. The overall message is, however, unlikely to be very different with Egypt-specific values for the impact of pollinator loss.
The total values of each of the Egyptian crops was obtained from the publication by the Economic Affairs Sector (2006) of the Ministry of Agriculture. This gives values either overall, or split by season or by land type (old or newly developed areas): we used the overall values.
The use of each crop, and therefore the impact of pollinator loss, varies. Thus some crops produce vegetative growth that is consumed: pollinators affect seed production for the next generation in those crops that are grown from seed each year. For crops that are grown vegetatively, the impact of pollinators is more long-term, but no less serious. For example, a standard fodder crop in Egypt is barseem (i.e. alfalfa), which can be cropped for six yearsbefore needing to be renewed from seed. However, alfalfa is well-known for its seed production being dependent on wild-bee pollinators because honeybees are especially poor: Brading et al.: semi-domesticated solitary bees (Megachile rotundata), on the other hand, do the job very well. Here we simply multiplied the value of the crop by the proportion of the yield that would be lost if pollinators were absent. Over the long term, we regard this as justifiable.
Results
The results (Table 1) are dominated (46%) by the impact of pollinator loss on alfalfa, predicted to cause annual losses of more than LE 6 billion (more than US $ 1 billion). It is true that even if seed production were reduced, this might not affect the production of the fodder itself.
However, in the long term, there would be a substantial impact. For valuable crops, such asmelons, pollinator loss would also have a huge economic impact, an annual loss of almost LE 1.9 billion (US $ 333 million).
Overall, according to these calculations, almost LE 13.5 billion (US $ 2.4 billion) would be lost every year. Since Egypt’s GDP in 2003 was LE 411 billion (UNESCO 2007), this represents about 3.3% of GDP.
Table 1: Monetary losses based on the average yield reduction that would be the consequence of loss of pollinators in Egypt for the crop production of 2004. There are approximately 5.7 LE to each $US. Pollination losses are taken from Klein et al (2007); crop values from Economic Affairs Sector (2006).
Commodity Latin name Pollination effect Pollination clover, alfalfa Trifolium spp, Medicago sativa Cotton Gossypium spp parts eaten 25 3131.3 782.8 Fodder (not alfalfa) various seeds 25 313.9 78.5 estimated % loss Lufa Luffa aegyptiaca seeds 65 109.8 71.4 probably 95% but some selfing occurs Sunflower Helianthus annuus parts eaten 25 96.5 24.1 Linseed, Flax, Straw Linum usitatissimum seeds 5 111.9 5.6 Safflower seed Carthamus tinctoria parts eaten 5 46.0 2.3 Kenaf Hibiscus cannabinus parts eaten 65 2.7 1.8 probably needs pollinators Egyptian lupin Lupinus albus seeds 5 11.4 0.6 mainly selfing Sugar Beet Beta vulgaris vulgaris seeds 0 357.6 0 Barley Hordeum spp independent 0 228.2 0 Rice Oryza spp independent 0 6678.6 0 Sugar Cane Saccharum officinarum independent 0 2191.1 0 Sorghum Sorghum spp independent 0 1001.8 0 Wheat Triticum spp independent 0 8903.9 0 Maize/corn/sweetcorn Zea mays independent 0 7361.4 0 Fruit crops Melon Cucumis melo parts eaten 95 1031.8 980.2 Mango Mangifera indica parts eaten 65 1323.3 860.1 Apple Malus 'domestica' parts eaten 65 930.7 605.0 Cantaloupe Cucumis melo parts eaten 95 547.7 520.3 Peach Prunus persica parts eaten 65 623.3 405.2 Watermelon Citrullus lanatus parts eaten 95 399.5 379.5 Apricot Prunus armeniaca parts eaten 65 169.0 109.8 Orange Citrus spp parts eaten 5 2166.6 108.3 Banana Musa spp breeding potential 5 1345.4 67.3 estimated % loss Fig Ficus carica parts eaten 25 258.4 64.6 Pear Pyrus communis parts eaten 65 91.3 59.4 Guava Psidium guajava parts eaten 25 235.9 59.0 Strawberry Fragaria spp parts eaten 25 134.6 33.6 Plum Prunus x domestica parts eaten 65 47.3 30.7 Tangerine, Mandarine Citrus spp parts eaten 5 595.1 29.8 Lemon, lime Citrus spp parts eaten 5 296.5 14.8 Prickly pears (Cactus) Opuntia parts eaten 25 48.9 12.2 Custard apple Annona spp parts eaten 95 10.8 10. Pomegranate Punica granatum parts eaten 25 39.5 9.9 Medlar (Loquat) Eriobotrya japonica parts eaten 65 2.1 1.4 Sour orange Citrus spp parts eaten 5 19.2 1.0 Kaki persimmon Diospyros kaki parts eaten 5 14.4 0.7 Grapefruit, Pomelo Citrus spp parts eaten 5 1.9 0.1 Olive Olea europaea independent 0 698.2 0 Date Phoenix dactylifera independent 0 1255.1 0 Grape Vitis vinifera independent 0 1912.5 0 Herb crops Rosemary Rosemarinus officinalis breeding potential 65 155.4 101.0 estimated % loss Marjoram Origanum majoranae breeding potential 65 58.6 38.1 estimated % loss Karkadeh Hibiscus sabdariffa parts eaten 65 51.4 33.4 probably needs pollinators Coriander Coriandrum sativum parts eaten 65 48.9 31.8 Basil Ocimum basilicum breeding potential 65 38.9 25.3 estimated % loss Wormwood Artemisia spp seeds 65 37.1 24.1 estimated % loss Fenugreek Trigonella foenumgraecum parts eaten 65 35.1 22.8 estimated % loss Parsley Petroselinum crispum breeding potential 65 27.9 18.2 estimated % loss Cumin Cuminum cyminum parts eaten 65 24.6 16.0 Sage Salvia spp seeds 65 22.9 14.9 estimated % loss Oregano Origanum vulgare breeding potential 65 22.9 14.9 estimated % loss Mint Mentha spp breeding potential 65 20.5 13.3 estimated % loss Fennel Foeniculum vulgare seeds 65 12.4 8.0 Anise Pimpinella anisum seeds 65 9.8 6.4 Dill Anethum graveolens seeds 25 24.1 6.0 estimated % loss Caraway Carum carvi parts eaten 25 17.5 4.4 Henna Lawsonia inermis parts eaten 65 2.8 1.8 estimated % loss Other aromatics various seeds 25 1.5 0.4 estimated % loss Coriander, green Coriandrum sativum seeds 65 0.4 0.3 Nut crops Almond Prunus dulcis parts eaten 65 201.2 130.8 Sesame seed Sesamum orientale parts eaten 25 157.4 39.4 Peanut, Groundnut Arachis hypogaea parts eaten 5 437.6 21.9 Pecan nut Carya illinoinensis independent 0 3.8 0 Vegetable crops Squash, courgette, pumpkin Cucurbita spp parts eaten 95 367.4 349.0 Cucumber Cucumis sativus parts eaten 65 348.7 226.6 Tomato Lycopersicon esculentum parts eaten 5 3797.0 189.8 Beans, Broad, dry Vicia faba parts eaten 25 757.9 189.5 Aubergine (eggplant) Solanum melongena parts eaten 25 398.4 99.6 Potato Solanum tuberosum breeding potential 5 1503.9 75.2 estimated % loss Beans, Broad, Green Vicia faba parts eaten 25 172.9 43.2 Molokhayia Corchorus olitorius seeds 95 43.6 41.4 grown from seed, and pollination required Okra Abelmoschus esculentus parts eaten 25 140.4 35.1 Onion Allium cepa seeds 5 580.7 29.0 estimated % loss Soybean Glycine max parts eaten 25 87.1 21.8 Carrot Daucus carota seeds 65 28.7 18.7 Snake Cucumber Cucumis melo parts eaten 65 27.3 17.7 Sweet peppers Capsicum annuum parts eaten 5 248.2 12.4 Artichoke Cynara scolymus breeding potential 25 35.7 8.9 estimated % loss Cabbage Brassica oleracea capitata seeds 5 174.3 8.7 estimated % loss Beans, dry Phaseolus spp parts eaten 5 142.6 7.1 Garlic Allium sativum breeding potential 5 139.4 7.0 estimated % loss Beans, green Phaseolus spp parts eaten 5 137.0 6.8 Sweet potato Ipomoea batatas breeding potential 5 88.2 4.4 estimated % loss Taro Colocasia esculenta breeding potential 5 72.2 3.6 vegetatively reproduced, but pollination by flies Radish Raphanus sativus parts eaten 65 4.9 3.2 annual, and mainly crosspollinated by insects Turnip Brassica rapa rapifera seeds 65 4.8 3.2 Rocket Eruca vesicaria sativa seeds 25 11.9 3.0 estimated % loss Beans, Kidney, Green Phaseolus spp parts eaten 5 49.8 2.5 Broccoli, Cauliflower Brassica oleracea botrytis seeds 5 43.2 2.2 estimated % loss Capsicum (chilli pepper) Capsicum annuum parts eaten 5 28.2 1.4 Onion seed Allium cepa seeds 5 23.1 1.2 estimated % loss Beans, Kidney, dry Phaseolus spp parts eaten 5 19.8 1.0 Egyptian leek Allium ampeloprasum var. kurrat seeds 5 9.1 0.5 estimated % loss Celery Apium graveolens seeds 5 2.3 0.1 estimated % loss Purslane (Rigla) Portulaca oleracea sativa seeds 25 0.3 0.1 estimated % loss Leek Allium ampeloprasum var. porrum seeds 5 0.4 0.02 estimated % loss Beetroot Beta vulgaris vulgaris independent 0 0.8 0 Chard Beta vulgaris vulgaris independent 0 6.8 0 Chick pea Cicer arietinum independent 0 29.5 0 Lettuce Lactuca sativa independent 0 42.8 0 Lentil Lens spp independent 0 8.8 0 Pea Pisum sativum independent 0 184.7 0 Spinach Spinachia oleracea independent 0 15.6 0 Total potential losses 13446.1
Discussion
Biological services, while essential for the whole planet’s survival and persistence, are often overlooked in a country’s economics. Without many of these services, invisibly working in the background, many economies would collapse. An early estimate for the value of pollination services was 0.4% of GDP for the whole world (Costanza et al. 1997); in managed pollination, a single solitary bee (Habropoda) can be worth $20 to Vaccinium pollination (Kevan & Phillips 2001), and the pollination services provided by nearby forest reserves for coffee plantations amounted to 7% of farm income (Ricketts et al. 2004). For comparison in the way in which we have calculated pollination services here, the value of insect ecosystem services to the USA was estimated by Losey & Vaughan (2006) at $58 billion (made up of dung burial 0.4, pollination 3.1, pest control of native herbivores 4.5, and ‘recreation’ [food for game, fish and wildlife] 50.0). Since the GDP of the USA in 2006 was $13 trillion, this represents only 0.45% of GDP, with pollination services accounting for only 0.02%.
In developing countries, pollinator services are almost certainly more significant in that a greater proportion of the human population is maintained by income provided by agriculture: Egypt is no exception. Crops that are fully dependent on pollinators, such as melons (including watermelons), onions and aubergines (eggplants), are some of the biggest contributors to the Egyptian agricultural market. With declining populations of pollinators, these crops will suffer a dramatic drop in production, and this will have a huge impact not only on the individual producers, but on the whole of Egypt’s economy.
The Nile Valley represents an environment with one of the world’s longest records of continuously habitation by man. Virtually all natural habitats have disappeared, and many insects must have been already lost before the advent of modern agriculture. 21st-century declines of pollinators on an already-narrowed group of pollinators are likely to be serious.
Egypt needs to implement strategies to prevent and reverse declines in pollinator populations.
Changing farming techniques (i.e. reducing intensification, conserving pollinator-friendly areas), and enforcing restrictions on pesticide use would go a long way to achieving this. While this sounds an expensive and counter-productive strategy, the potential consequences of not implementing such a change could be far more costly to Egypt’s development.
Acknowledgements: We thank Dr Moustafa Fouda for facilities in the Nature Conservation Sector of the Egyptian Environmental Affairs Agency, and for reviewing the manuscript.

El-Gabbas A. & Dormann C. F.: Improved species-occurrence predictions in data-poor regions: using large-scale data and bias correction with down-weighted Poisson regression and Maxent. Ecography DOI: 10.1111/ecog.03149. Appendix 1: Supplementary figures and tables Table A1: The estimated optimum combination of Maxent’s feature classes (FC) and Regularization Multiplier (RM) for each species and bias model type. Combinations with highest mean testing-AUC on 5-folds spatial-block cross-validation were selected. All the analyses were performed using modified code from the ENMeval package in R (Muscarella et al., 2014). See main text for more information. Species Environment-only Accessibility Effort a Feature classes (FC): L linear; Q quadratic; H hinge; P product; and T threshold. For more details see Elith et al. 2011. b RM is a global multiplier to all individual regularization values; for more details see the main text and Muscarella et al. (2014). Table A2: Estimated best values to run the elastic-net models (DWPR approach) for different species and bias models. For details, see main text. Species Environment-only Accessibility Efforts LASSO RIDGE Fig. A1: Map of all non-marine GBIF-records. There are clear signs of sampling bias inherited in the GBIF database, disallowing the direct use of these records to run SDMs, unless action is taken to correct for sampling bias. GBIF-records are biased mainly towards Western Europe, with sparse locations across Africa, Asia, and Eastern Europe. Areas shown in white represent pixels without any records (the majority of North and central Africa, Arabia, and eastern Europe). Clearly errorneous records were excluded before plotting this map.Fig. A2: Left: Map of all GBIF bat records (accessed date: April 2015) per a grid of 20 × 20 km2.Right: All records available of the study species, collected either from the literature (blue points) or GBIF (red points). Both maps show large collection gaps. Fig. A3: Per-species study area: An example of how the ‘per-species’ study area was determined (here for Asellia tridens). The dashed red line indicates the species extent of occurrence (minimum convex polygon); the solid red line shows the buffer used (1000 km), and the dashed green line indicates the study area used. Fig. A4: Variables selection: Pearson correlation coefficient between each pair of the final covariates; see Appendix 3 for more information. Fig A5: Per-species number of backgrounds: Objectively determining the adequate number of background ‘quadrature’ points to be used in GLM and elastic-net models (see main text and Renner et al., 2015). DWPR-GLM models (for each species, here: Asellia tridens) were repeated 25 times (each time, with a different random sample from the background) with progressively increase the number of background points (from 5000 to 500,000); and for each model, the log-likelihood value is calculated and plotted. The number of backgrounds at which the likelihood converges (i.e. no much benefit of using a higher number of backgrounds) is, visually, selected to run the final models. Here, no much improvement in the log-likelihood beyond 300,000 background points (out of 557,800 total available background points) and so we used it for this species. Fig. A6: Determining the best value to run the elastic-net model for Nycticeinops schlieffeni (effort model). Eleven values ranging from zero (ridge) to one (lasso), with an increment of 0.1, were used to run 11 models on cross-validation. For each value of , glmnet ran many models on a range of -values (estimated from the data; x-axis) and measure the mean cross-validation error (cvm) (here: Poisson deviance; y-axis). The value with the lowest cvm is used further while predicting (vertical dotted lines, with colours correspond to the value used). Here, an value of 0.6 is shown to have the lowest cvm and so used further. For more details, see Table A2 and the main text. Fig. A7: Kendall’s correlation between (per-species and modelling algorithm) mean AUC evaluated on spatial-block cross-validation and independent evaluation in Egypt. Each point represents the mean of 100 AUC values, with different symbols for each species. Colours represent different modelling algorithm used. ‘M’ indicates the overall mean for each modelling algorithm. Panel (a) shows the correlation for the environment-only models (no bias correction). Panels (b-c) represent the accessibility model, without and with correcting for sampling bias, respectively. Panels (d-e) represent the effort model, without and with correcting for sampling bias, respectively. Kendall’s correlation coefficients (and their p-values) are reported in each block. Grey solid line represents the equality line. AUC performed on cross-validation were, on average, higher than AUC in Egypt. However, there is moderate correlation between bias-free evaluations on either scale. Fig. A8: Equivalent results to Fig. A7, for TSS evaluation. Fig A9: Boxplots of the raw cross-validation evaluation (100 values) without (modelling evaluation; a and c) and with (bias-free evaluation; b and d) correction for sampling bias; either using AUC (a - b) or TSS (c - d); see Appendix 6 for more information. Horizontal panels show results for different modelling algorithm, while the vertical panels show different bias models. Panels a & c show results for the environment-only model (without any bias manipulation) and bias-accounting models (accessibility and effort models) without correcting for sampling bias. However, panels b & d compare results for the environment-only model and bias-accounting models (accessibility and effort models) after correcting for sampling bias (see main text; for effort models, plots show evaluations either conditioning the predictions on a value of zero or the maximum relative sampling effort of training presences). Species are in ascending order according to their number of occupied pixels at cross-validation scale (with numbers represent the species; see Table 1 for full species names). Fig. A10: Similar to Fig. A9, but showing results of independent evaluation in Egypt. Evaluation metrics are calculated in Egypt using mean predictions (of 5-folds cross-validation) in Egypt along with entirely independent species records from Egypt (not used to run any of the models). Species are in ascending order according to their number of occupied pixels in Egypt (with numbers represent the species; see Table 1 for full species names). Fig. A11: Species mean TSS calculated either on cross-validation (a - b) or in Egypt (c - d), either without (modelling evaluation; a & c) or with (bias-free evaluation; b & d) sampling bias correction (for details, see Appendix 6). Each species is represented by different symbols (similar to those shown in Fig. A7; numbers represent species used, see Table 1). Red lines indicate the overall mean TSS at each modelling algorithm and bias models applied. Fig. A12: Kendall’s correlation of the per-species mean TSS between different pairs of modelling algorithms (same as Fig. 2 in the main text, but for TSS). Each species is represented by different symbol (similar to those of Fig. A7) with different colours for different bias models applied (using predictions of environment-only model and bias-free prediction of accessibility and effort model). ‘M’ indicates the overall mean evaluation. First row represents mean spatial-block cross-validation evaluation, while the second is for independent evaluations in Egypt. Fig. A13: Values of bias variables at presences and available locations (backgrounds), at the per-species study area (A) or at local scale ‘Egypt’ (B). Rows correspond to different species, with numbers indicates the species (see Table 1), and columns represent bias variables of the accessibility model (distances to roads, cities, and protected areas) and effort model (relative bats’ sampling intensity). For each species, values at presence locations are indicated with black, and values at pixels unoccupied by the species are shown in grey. Most of the species are recorded from closer to roads and cities (and to some extent, the protected areas), and unexpectedly at low to moderate sampling efforts. Fig. A14: AUC scores calculated the default way (using all available testing backgrounds – y-axis) versus using a fixed ratio between testing presences and backgrounds (test-data prevalence = 1:20 – x-axis). (A) shows the raw outputs of 5-folds cross-validation. (B) shows the per-species mean AUC on cross-validation. Different colours represent the three bias models applied (environment-only, accessibility, and efforts). Fig. A15: The predicted distribution of Otonycteris hemprichii (mean of 5-folds cross-validation), using different modelling algorithms (rows) and bias models (columns). Maps were scaled between zero and one, as different modelling algorithms do not have the same scale, with blue colour indicates higher predicted relative intensity. (A) shows cropped predicted distribution of Otonycteris hemprichii (the same as Fig. 5) to Egypt. Grey points (in the top left panel) represents available records used for indepenent evaluation presences in Egypt. In Fig. 5 of the main text, extreme values (> 0.9995 quantile of predicted values) were replaced with their next smaller value to improve map visualization, as GLM and elastic-net models sometimes yielded extreme values. Maps in (B) are equivalent to maps in Fig. 5, without any extreme values manipulation (using the linear scale), demonstrating the difficulty of visualizing the predicted patterns in the existence of extreme values; see main text for more details. Fig. A16: The reported and predicted distribution of Rhinolophus mehelyi (Maxent). The 2nd to the 4th panels show predictions from environment-only model and bias-free prediction of accessibility and effort models, respectively. Arrows represent regions which gain higher predicted values after correcting for sampling bias (such as central Turkey or the Algerian Atlas). Predictions at these areas are of higher uncertainty (lower congurrence) and field validation is probably required to confirm the species existence. The lighter the colour, the more suitable the location. Appendix 3: Variables selection List of initial covariates investigated for correlation/multi-collinearity and transformation applied. Final list used to run the model is marked in grey. Variable Transfo- rmation VIF 10 | cor | 0.7 GVIF 3 Alt a Altitude yes yes yes Bio1 Annual mean temperature ^2 Bio2 Mean diurnal range (mean of monthly (max temp - min temp)) yes yes yes Bio3 Isothermality (Bio2/Bio7) (* 100) Bio4 Temperature seasonality (standard deviation *100) yes yes yes Bio5 Maximum temperature of warmest month Bio6 Minimum temperature of coldest month Bio7 Temperature annual range (Bio5-Bio6) Bio8 Mean temperature of wettest quarter ^2 yes yes yes Bio9 Mean temperature of driest quarter ^2 yes yes yes Bio10 Mean temperature of warmest quarter ^2 Bio11 Mean temperature of coldest quarter yes yes Bio12 Annual precipitation sqrt Bio13 Precipitation of wettest month yes Bio14 Precipitation of driest month yes yes yes Bio15 Precipitation seasonality (coefficient of variation) yes yes yes Bio16 Precipitation of wettest quarter sqrt Bio17 Precipitation of driest quarter Bio18 Precipitation of warmest quarter yes yes Bio19 Precipitation of coldest quarter log yes yes yes PET b Potential evapotranspiration AI b Aridity index sqrt AET c Actual evapotranspiration sqrt SWB c Soil-water balance sqrt NDVI_Max d Maximum NDVI (Normalized Difference Vegetation Index) NDVI_Min d Minimum NDVI yes yes yes NDVI_Mean d Mean NDVI NDVI_Range d Range NDVI (maximum – minimum) NDVI_SD d Standard deviation of NDVI yes yes yes EVI_Max d Maximum EVI (Enhanced Vegetation Index) EVI_Min d Minimum EVI yes EVI_Mean d Mean EVI yes EVI_Range d Range EVI (maximum - minimum) yes EVI_SD d Standard deviation of EVI a WorldClim : WorldClim provides a global dataset for the elevation and 19 bio-climatic variables interpolated from global monthly temperature and precipitation recordings from weather stations (Hijmans et al. 2005). Tiles for the overall study area were downloaded at 30 arc-seconds resolution (~1 km near the equator; using the raster R_package), then projected to Mollweide equal-area projection at 5 × 5 km2 resolution. The high resolution of 30 arc-seconds was preferred over the 2.5 arc minutes (~5 km near the equator) in order not to lose much information while re-projecting. Instead of interpolation (assigning a value for the new pixel equals to the mean of the nearest three points at the original projection), a different approach was employed. First, a template grid covering the study area at the equal-area projection (5 × 5 km2) was prepared. Then, pixels of the original variables (at 30 arc-seconds resolution) were converted into centroid points (at the original projection: WGS-1984), then projected into Mollweide projection, then these points were rasterized using the mean value (or other relevant function for some variables; e.g. maximum for Bio5) of the points that spatially fall within each cell of the template grid. The same approach was used to prepare all other covariates used in this study. b Global Potential Evapotranspiration (Global-PET) & Global Aridity (Zomer et al. 2007, 2008) : available at the global scale at a resolution of 1 km, and were prepared based on models that use data from WorldClim as inputs. Global Actual Evapotranspiration (Global-AET) & Global Soil-Water-Balance (Global-SWB) (Trabucco & Zomer 2010) : available at the global scale at a resolution of 1 km, and are prepared based on WorldClim and Global-PET database as primary input. d NDVI (Normalized Difference Vegetation Index) & EVI (Enhanced Vegetation Index) (Didan, 2015) : map tiles for the whole study area were downloaded (and merged) for the period from 18/2/2000 to 22/3/2015 each 16 days (MODIS product: MOD13A2 – resolution: 1 km) using the MODIS R package (Mattiuzzi, 2014). Summary maps (maximum, minimum, mean, range, and standard deviation) across the whole period (for NDVI & EVI) were produced at the original MODIS scale/projection using a python code in ‘ArcGIS’ as this was memory intensive to perform in raster package of R. Each of the summary layers were then projected into the equal-area projection at 5 × 5 km2 resolution. Possible covariate transformations were investigated, trying to keep some degree of uniform distributions across the range of each covariate (following: Dormann 2011). For some covariates, no transformation was effective, and they were kept untransformed. • Variable selection: Multicollinearity amongst the potential covariates )covering the total study area( was assessed, maintaining highest GVIF (Generalized Variance Inflation Factor) 3 (Zuur et al. 2009) [which is equivalent, in our case, to maintaining a highest absolute correlation coefficient between each pair of covariates 0.6 and highest VIF (Variance Inflation Factor) 10]. This was also checked for each species’ study area. Aridity-related covariates (actual/potential evapotranspiration, aridity-index, and soil-water-balance) were all excluded as they show high collinearity with other WorldClim covariates; unsurprisingly as they were derived from models that use WorldClim data as input. Out of tensummary vegetation-related covariates, two NDVI covariates were used further: the cumulative minimum and the standard deviation NDVI. They reflect indices of minimum biomass and variability of vegetation-cover across the study area, respectively. Appendix 4: Per-species spatial-block cross-validation For each species, a different spatial-block cross-validation structure was used to run SDMs. Pixels across the study area were aggregated into larger blocks (each of 20 × 20 cells = 100 × 100 km2), and blocks were distributed into cross-validation folds. Presences and backgrounds in each block were used together either as training or testing (Fithian et al. 2015). The blocks were not distributed into cross-validation folds randomly, as this would have yielded highly unbalanced numbers of presence-locations across spatial folds. Instead, we balanced the number of presence locations at different folds by calculating their numbers at each block; the top 5 blocks were then, sequentially, randomly assigned to five folds (to avoid that the first fold always has the highest number of presences). To avoid high variability in the environmental space between folds and minimise much extrapolation while predicting to the left-out-fold, blocks without any species records were distributed into folds depending on a mean index of their similarity to the overall study area: using the Multivariate Environmental Similarity Surfaces ‘MESS’. MESS was proposed by Elith et al. (2010) to identify areas of novel climates by comparing future (or past) climate to those used to run the models, as the interpretation of future projections at these locations should be considered with caution. We calculated the MESS-score for each pixel (and a mean value for each block) quite differently from how it was proposed in the first place: we measured the environmental similarity of each pixel to all available pixels in the study area (using the ‘dismo’ package). The original application of MESS gives most dissimilar locations more negative values (Elith et al. 2010); however here we expect no negative values as there is no climates novelty. Blocks without any presence locations are sorted in descending order (depending on their average MESS value), and then the highest five blocks were distributed randomly into different folds. This was repeated for all blocks until all blocks were distributed to different folds. An example of how spatial-block cross-validation applied for Asellia tridens. Different colours represent how different blocks are distributed into cross-validation folds. Blue points represent the species known distribution from outside Egypt. For each per species’ study area, a different blocking structure was prepared, to balance the number of presence locations and environmental variability at different cross-validation model runs. Appendix 5: Sampling-effort model Data for closely related species (bats/non-marine mammals) were downloaded from GBIF (The Global Biodiversity Information Facility – April 2015) and used (along with available focal species records: Tables 1 & A1) to model the relative sampling intensity (surveying effort) of bats/mammals. The prediction map of this model was used as bias covariate in the ‘Effort models’ (see main text, Fig. 1, bottom right). Records were assessed before usage (records with missing coordinates or with clear errors were excluded: e.g. missed latitude or longitude / equal latitude and longitude / records in the sea or the ocean / potentially swapped latitude and longitude), resulting in a total of 435,458 bat records / 2,039,158 non-marine mammal records (maps below). GBIF records (for both bats and mammals) show high bias towards the Western Europe compared to any other area in the study area, followed by scattered locations elsewhere (mostly close to the main cities, water bodies, populated areas, or seemingly a result of a mammals atlas mapping activity in southwest Africa). The majority of Africa and eastern Europe to western Asia is extremely under-represented in GBIF, with huge gaps in North Africa and Arabia. Maps show the number of GBIF records of non-marine mammals (left) and bats (right) per a grid of 20 × 20 km2 (accessed date: April 2015). Initial trials for modelling the sampling effort were done using DWPR-GLM (Down-Weighted Poisson Regression) and Maxent (both using PPM-like approach; see main text); without much difference in the resulted prediction pattern (although Maxent model ‘using default feature classes and regularization multiplier’ shows, visually, an over-fitted spatial pattern). DWPR-GLM model was chosen to run at 5×5 km2 equal-area projection: all the bats/mammals presence locations as the response (without duplicates removal) and using different covariate set than the environmental variables used to run the species SDMs (Merow et al. 2016) Terrain roughness (represented here as a per-pixel standard deviation of altitude): Altitude maps from WorldClim was downloaded at of 30 arc-seconds resolution [~1 km near the equator], then projected as points in Mollwide equal-area projection. The points were then converted to a raster (rasterized), representing the standard deviation of the elevation values located at any target pixel at the resolution of 5 × 5 km2; Distance to main cities : the Euclidean distance between the centroid of each pixel and the nearest human settlement; Distance to main roads : the Euclidean distance between the centroid of each pixel and the nearest road (Source: Global Roads Open Access Data Set ‘gROADS’); Population count : global population count in 2000; Protection status : a binary variable indicating the protection status of each pixel [source: The world Database of Protected Areas:]. The overall pattern of the predicted sampling effort for both bats and mammals were roughly similar (maps below), so the prediction from the bats sampling effort model was used further (as bias covariate in the effort model). It represents the overall relative abundance of bat species sightings per unit area (based only on non-climatic covariates). The predicted sampling intensity of sightings (sampling effort) for all bats (left) and non-marine mammals (right) using the DWPR-GLM (Down-Weighted Poisson Regression) model. Appendix 6: Mixed effect model analysis of evaluation Mixed Models formula: Evaluation metric ~ Modelling algorithms * Bias models + total number of training presences (training data) + number of pixels in the study area (total area) + (1 | species) Evaluation metric continuous variable for either AUC or TSS either using mean cross-validation evaluation (AUCcv/TSScv) or independent evaluation in Egypt (using mean prediction of 5-folds cross-validation to Egypt and independent testing data from Egypt) Modelling algorithms: categorical variable indicating the modelling algorithm used (elastic net / GLM / Maxent) Bias models: 1. Modelling evaluation: evaluation of environment-only model (base line; no bias manipulation) is compared with those of bias-accounting models (accessibility and effort), without correction for sampling bias [without fixing bias variables at any values]. 2. Bias correction evaluation: evaluation of environment-only model (base line; no bias manipulation) is compared with those of bias accounting models (accessibility and effort), after correction for sampling bias [accessibility bias variables are set to zero, while the effort bias variable is set to either zero or the maximum estimated effort of the target species’ presence records; see below]. Species: random factor; categorical variable for the species. Modelling Evaluation – Without Bias correction Variance explained (sum of squares) AUC TSS Mean cross-validation Egypt evaluation Mean cross-validation Egypt evaluation Bias models Model techniques Bias models * Model techniques Total number of training presences Study area Parameter estimates (Mixed models summary) Bias correction evaluation – using bias-free predictionsAfter correcting for sampling bias, initial trials show that using either value for fixing the effort bias covariate produce very similar evaluations (however, Maxent models show quite lower AUCs using the maximum effort value of training presences than for fixing at zero; Fig. 3), so we limited further analyses to effortMax.When predicting with fixed values for the bias-covariates, modelling algorithms were most important for explaining variability of cross-validations, followed by sampling-bias models and then their interaction (for evaluation in Egypt highest variability was explained by the sampling-bias models, followed by the modelling algorithm, and then their interaction). Again, the number of training presences and study area were much less important (and the sign of their effects resembles those of modelling evaluation). The accessibility model had lower AUC compared to the two other sampling-bias models in all comparisons (Fig. 3). On cross-validation, bias-accounting models had lower AUC-values compared to the environment-only model (very little difference between environment-only and effort model for GLM; Fig. 3a). For evaluations in Egypt, environment-only and effort model were hardly different (effort model had lower AUC for Maxent ; Fig. 3b).Variance explained (sum of squares) AUC TSS Mean cross-validation Egypt evaluation Mean cross-validation Egypt evaluation Bias models Model techniques Total number of training presences Study area Parameter estimates (Mixed models summary) AUC TSS Mean cross-validation Egypt evaluation Mean cross-validation Egypt evaluation(Intercept) [env-only/ elastic net] Estimated differences in the least square means AUC TSS Mean cross-validation Egypt evaluation Mean cross-validation Egypt evaluation Modelling techniques Supplementary references Didan K (2015): MOD13A2 MODIS/Terra Vegetation Indices 16-Day L3 Global 1km SIN Grid V006. NASA EOSDIS Land Processes DAAC. Dormann CF (2011): Modelling Species’ Distributions. In Fred J, Hauke R, Broder B (Eds.): Modelling Complex Ecological Dynamics. Berlin, Heidelberg: Springer Berlin Heidelberg, pp. 179–196. Elith J, Kearney M, Phillips S (2010): The art of modelling range-shifting species. Methods in Ecology and Evolution 1: 330–342. Fithian, W., Elith, J., Hastie, T., Keith, D.A. & O'Hara, R.B. (2015). Bias correction in species distribution models: pooling survey and collection data for multiple species. Methods in Ecology and Evolution, 6, 424–438. Friedman J, Hastie T, Tibshirani R (2010): Regularization Paths for Generalized Linear Models via Coordinate Descent. Journal of Statistical Software, 33(1), 1-22. URL: http://www.jstatsoft.org/v33/i01/. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A (2005): Very high resolution interpolated climate surfaces for global land areas. In Int. J. Climatol. 25 (15), pp. 1965–1978. Merow C, Allen JM, Aiello-Lammens M & Silander JA (2016): Improving niche and range estimates with Maxent and point process models by integrating spatially explicit information. Global Ecology and Biogeography, 25, 1022–1036. Mattiuzzi M (2014): MODIS: MODIS acquisition and processing. R package version 0.10-17/r484. (https://R-Forge.R-project.org/projects/modis/) Renner IW, Elith J, Baddeley A, Fithian W, Hastie T, Phillips SJ., et al. (2015): Point process models for presence-only analysis. Methods in Ecology and Evolution 6: 366–379. Trabucco A, Zomer RJ (2010): Global Soil Water Balance Geospatial Database. CGIAR Consortium for Spatial Information. Published online, available from the CGIAR-CSI GeoPortal at: http://www.cgiar-csi.org. Zomer RJ, Bossio DA, Trabucco A, Yuanjie L, Gupta DC, Singh VP (2007): Trees and Water: Smallholder Agroforestry on Irrigated Lands in Northern India. Colombo, Sri Lanka: International Water Management Institute. pp 45. (IWMI Research Report 122). Zomer RJ, Trabucco A, Bossio DA, van Straaten O, Verchot LV (2008): Climate Change Mitigation: A Spatial Analysis of Global Land Suitability for Clean Development Mechanism Afforestation and Reforestation. Agric. Ecosystems and Envir. 126: 67-80. Zuur AF, Ieno EN, Walker NJ, Saveliev AA, Smith GM (2009): Mixed Effects Models and Extensions in Ecology with R. New York, NY: Springer New York.

El-Gabbas A. & Dormann C. F.: Improved species-occurrence predictions in data-poor regions: using large-scale data and bias correction with down-weighted Poisson regression and Maxent. Ecography DOI: 10.1111/ecog.03149. Ahmed El-Gabbas & Carsten F. Dormann Department of Biometry and Environmental System Analysis, University of Freiburg, D-79106 Freiburg, Germany bat, Maxent, point-process model, presence-only data, species distribution modelling, sampling bias Running-title: Bias correction for sparse presence-only data Abstract Species distribution modelling (SDM) has become an essential method in ecology and conservation. In the absence of survey data, the majority of SDMs are calibrated with opportunistic presence-only data, incurring substantial sampling bias. We address the challenge of correcting for sampling bias in the data-sparse situations. We modelled the relative intensity of bat records in their entire range using three modelling algorithms under the point-process modelling framework (GLMs with subset selection, GLMs fitted with an elastic net penalty, and Maxent). To correct for sampling bias, we applied model-based bias correction by incorporating spatial information on site accessibility or sampling efforts. We evaluated the effect of bias correction on the models’ predictive performance (AUC and TSS), calculated on spatial-block cross-validation and a holdout data set. When evaluated with independent, but also sampling-biased test data, correction for sampling bias led to improved predictions. The predictive performance of the three modelling algorithms was very similar. Elastic-net models have intermediate performance, with slight advantage for GLMs on cross-validation and Maxent on hold-out evaluation. Model-based bias correction is very useful in data-sparse situations, where detailed data are not available to apply other bias correction methods. However, bias correction success depends on how well the selected bias variables describe the sources of bias. In this study, accessibility covariates described bias in our data better than effort covariates, and their use led to larger changes in predictive performance. Objectively evaluating bias correction requires bias-free presence-absence test data, and without them the real improvement for describing a species’ environmental niche cannot be assessed. Introduction Species distribution data often come in the form of presence-only, with information on where species have been recorded, but no reliable information on where they have not, or where people have looked (Pearce and Boyce 2006). Museums, herbaria, personal collections, published literature, and citizen records are valuable sources for presence-only data (Pearce and Boyce 2006, Newbold 2010), especially in developing countries where there is a lack of systematic nation-wide surveys. Some initiatives make such species sightings freely available: GBIF (the Global Biodiversity Information Facility – see: http://www.gbif.org) collates global biodiversity data from different sources. One fundamental problem is that data are often incidental with no information on the sampling efforts and survey method used (Pearce and Boyce 2006). They are biased taxonomically (towards larger, easy to detect, or more charismatic species groups; Newbold 2010), environmentally (less collection effort in areas with harsh environments), temporally (more in summer than winter), and spatially (near populated places, roads, research institutes and protected areas; Phillips et al. 2009, Newbold 2010, Stolar and Nielsen 2015). For example, GBIF-data show huge differences in data contribution among countries (Supplementary material Appendix 1, Fig. A1; on average, more from well-financed than from species-rich countries). Spatial bias is a particular concern for statistical analysis when it leads to environmental bias (Phillips et al. 2009), e.g. when large parts of the environmental space remain unsampled (Merow et al. 2014). Species distribution models (SDMs) relate species occurrences to the environment to estimate habitat preference, and predict potential distribution and responses to climate change (Phillips and Dudík 2008, Elith et al. 2011). Statistical analyses of presence-only data describe the environment at record locations relative to the background environment, making them more susceptible to sampling bias than presence-absence data from dedicated surveys (Phillips et al. 2009, Fithian et al. 2015). However, such targeted survey-data are rare (Pearce and Boyce 2006), especially in developing countries, explaining why the majority of SDM applications uses presence-only data. Presence-only data may produce sound models if they are efficiently corrected for sampling bias (Elith et al. 2011). Point-process models (PPM) have recently emerged as the most appropriate technique for presence-only data (Renner et al. 2015). PPMs do not use background points as pseudo-absences (as in the naïve logistic regression; Fithian and Hastie 2013), but rather as quadrature points for estimating the spatial integral of the likelihood function, and hence require careful tuning of its number (for details see, e.g., Warton and Aarts 2013). The response variable of PPMs is the density of species records per unit area (also called ‘intensity’), which should be proportional to the probability of occurrence (which can not be estimated empirically using presence-only data without additional information; see: Fithian and Hastie 2013, Renner et al. 2015, Phillips et al. 2017). It is mathematically equivalent to methods already commonly used in ecology, e.g. Maxent and some implementations of the generalised linear modelling framework, but differently efficient (Renner and Warton 2013, for details, see: Renner et al. 2015). Sampling bias has been addressed by ‘spatial filtering’ of presence locations (keeping only a limited number of records within a certain distance) to dilute the effect of uneven sampling effort across the study area (e.g. Anderson and Raza 2010, Boria et al. 2014). Alternatively, others effectively use location of records from related species as background points to have background points with the same bias as the species records (Elith and Leathwick 2007, 'target-group background': Phillips et al. 2009, see also: Ponder et al. 2001, and 'weighted target group' presented in Anderson 2003). Neither approach is applicable when there are only few data (a typical case in developing countries). For example, to apply the target-group background approach to model the distribution of a bat species in North Africa, all GBIF bat species records (Supplementary material Appendix 1, Fig. A2) are clearly not enough to serve as representative background points in this large study area. Similar to the target-group background, if the pattern of sampling bias is known a priori, it can be used as prior weight for sampling the background proportionally to the sampling effort (e.g. 'bias file' in Maxent; Phillips and Dudík 2008, Warren et al. 2014), so that both presences and background samples have the same bias (see also: Stolar and Nielsen 2015). As a third strategy, sampling bias can be addressed also by modelling the distribution of the focal species as a function of two additive covariate sets: the environmental covariates and other covariate(s) describing potential sources of sampling bias, hereafter ‘bias covariates’ (model-based bias correction; Warton et al. 2013). For unbiased predictions, the bias covariates are set to a common level of bias, say 0, at all locations; however, sometimes it is difficult to settle on a meaningful adjustment level (Warton et al. 2013). The aim of this paper is to address the problem of sampling bias in data-sparse situations and how to correct for it in presence-only SDMs. We apply model-based bias correction, comparing two different sets of bias covariates to model the distribution of Egyptian bat species in areas of their known global distributions. Bias-models for each of three modelling techniques within the PPM framework are calibrated with information on either accessibility or sampling effort and compared to an ‘environment-only model’. To maintain a reasonable degree of independence between training and testing datasets, we evaluated the models using a) entirely independent presence-only data (not used to fit any of the models); and b) spatial-block cross-validation.   Materials and methods Species and study area We are interested in understanding the environmental preferences of Egyptian bats as an example representing the sampling-bias issue in data-sparse situations. We collected records for the entire range of each species (Supplementary material Appendix 1, Fig. A2, latitude: −35° to +56° and longitude: −20° to +80°) from the literature and GBIF (see Supplementary material Appendix 2 for list of literature sources). Although GBIF provides a valuable source of data, such opportunistically compiled data bases inevitably contain misidentified or incorrectly georeferenced records (Gaiji et al. 2013), and thus require careful revision before use. Relevant records from the GBIF database were assessed (October 2014) and merged with the available literature records. Bat records from Egypt were mostly taken from the expert-revised BioMAP (Biodiversity Monitoring and Assessment Project) database (Basuony et al. 2010). Only species with enough presence locations that are located in at least 5 larger spatial blocks were included in this study (allowing model evaluation on spatial-block cross-validation; see below). This was fulfilled by 21 species, each occupying > 20 unique cells at resolution of 5×5 km2 (Table 1). The geographical range of records was assessed (based on the literature and IUCN), and spatial outliers were excluded. The coverage of all available records shows obvious signs of spatial bias towards Western Europe and only sparse sampling in Africa and Western Asia (Supplementary material Appendix 1, Fig. A2). Presence locations were purposefully split into training and testing data, using only records from outside Egypt’s boundaries for training, thereby keeping the Egyptian presence data as entirely independent evaluation data. For sampling of background points, however, Egypt was not excluded from the study area, and hence background points can be sampled from Egypt as well (see Fig. 1 for a flowchart of the methods applied). The determination of the study area is critical, especially for presence-only models (Pearce and Boyce 2006). We decided against a single fixed large study area that covers presence locations of all 21 species to keep only areas of potential accessibility to the bats (Barve et al. 2011) and avoid inflating the discrimination ability of the model (e.g. higher AUC; Barve et al. 2011). Instead, for each species, the study area was determined based on the geographical extent of the records: a rectangular bounding box containing a 1000 km buffer around the species extent of occurrence (see Supplementary material Appendix 1, Fig. A3 for an example). We used a buffer of 1000 km, as we found it suitable for the study species ‘bats’, which have, on average, high dispersal ability and large home range. Potential covariates (and species presences) were projected into Mollweide equal-area projection and rasterised to a resolution of 5 × 5 km2. We assessed potential environmental variables for multi-collinearity, and assembled a final list of covariates with a maximum Generalized Variance Inflation Factor value less than 3 (for details see Supplementary material Appendix 3 & Supplementary material Appendix 1, Fig. A4). Block cross-validation Cross-validation is commonly used for evaluating model performance when no independent data are available. Random splitting does not guarantee spatial independence due to spatial autocorrelation (both training and testing data will be spatially adjacent), and so may overestimate model performance (Bahn and McGill 2013, Radosavljevic and Anderson 2014). To maintain independence between folds and improve transferability of models, a spatial form of cross-validation was used by splitting the study area into coarse checkerboard blocks, and then distribute blocks randomly into folds (Fithian et al. 2015). We performed 5-fold spatial-block cross-validation (hereafter: cross-validation). The larger the block size, the higher is the need for more data to be able to run models on spatial-block cross-validation. We used blocks of 100 × 100 km2 (20 × 20 cells), which is not very strong, as we find this more appropriate for the available data. Presence and background locations within each block were used together for model training or testing: for the three modelling algorithms, potential background locations of the left-out cross-validation fold were not used (masked) during model calibration and were used exclusively for evaluation (similar to 'masked geographically structured approach' of Radosavljevic and Anderson 2014, see also Fig. 7 in Fithian et al. 2015). For each species, we used a different blocking structure, balancing the number of available presences between folds and avoiding extrapolation in environmental space (for more details see Supplementary material Appendix 4). Sampling-bias models We compared models without bias correction (environment-only model, our reference) with two methods of addressing bias. Firstly, we considered human accessibility as the main source of sampling bias and ran SDMs with additional bias covariates describing distances to nearest cities, roads and protected areas (the ‘accessibility model’). For bias-free predictions, we set all distances to zero; thus, a bias-free prediction can be interpreted as the relative intensity of the target species records if all locations across the study area had perfect accessibility (Warton et al. 2013). Second, assuming relative effort is the main source of sampling bias, we incorporated a single bias covariate describing the relative intensity of sightings of all bat species (the ‘effort model’). This sampling-effort bias covariate was actually the prediction from a (different) all-bats-model, which predicted the number of all records as a function of non-environmental variables (terrain roughness, distance to cities, distance to main roads, human population density, protected area – details in Supplementary material Appendix 5 and Fig. 1 bottom-right). For bias-adjusted predictions, the sampling-effort covariate needed to be set to a common level, for which there is no obvious choice. We therefore adjusted it to either of two values of sampling effort: the maximum fitted value at the target species' training presence locations, and at zero. In the first case, the bias-adjusted prediction can be interpreted as the relative intensity of species reporting if all the study area receives the sampling effort of the best presence location of the target species. In the second case, the bias-adjusted prediction is independent of effort, and the effort covariate is used only to correct the coefficients of the environmental variables. Although the latter value seems simpler, the predicted values lack intuitive interpretation. As the results were very similar, we here only use the maximum fitted value at training presences. Modelling algorithms For each sampling bias model (environment-only, accessibility and effort), we employed three modelling algorithms under the PPM framework with cross-validation and adjustment for model complexity: GLM, elastic net, and Maxent. Both GLM and elastic net model the number of records as a Poisson regression, with elastic net including a mixture of lasso and ridge regularisation (shrinkage; L1- and L2-regularisation, respectively; Friedman et al. 2010). Maxent (v3.3.3k; Phillips and Dudík 2008) is a machine-learning algorithm for presence-only data, effectively and mathematically akin to a Poisson GLM, but with different functional forms for the predictors; it also applies an ad-hoc form of lasso regularisation (Renner and Warton 2013). We used GLM and elastic net to implement a “down-weighted Poisson regression” (DWPR, following Renner et al. 2015). This approach uses weights to make the number of presences a negligible proportion of all data and scales the data to the actual area. As a consequence, DWPR estimates model parameters including the intercept. A small weight (10-6) was assigned to presence locations, while background points were given a higher weight equal to the area of the study region divided by the number of background points used. To estimate the appropriate number of background points (for GLM and elastic net), 25 repeated series of DWPR-GLMs were run, each one progressively increasing the number of randomly sampled background points. We used the number of background points at which the log-likelihood asymptoted (Supplementary material Appendix 1, Fig. A5; Renner et al. 2015). For both GLM and elastic-net models, 1) we included linear, quadratic and 1st-order interactions between variables; 2) no interactions were allowed between the environmental and bias covariates (Warton et al. 2013); 3) we standardised all covariates to a mean of zero and standard deviation of one. For GLMs, a random sample of background points (number estimated from the asymptoted log-likelihood curve) was used to run an initial model. The initial model was then simplified using AIC-informed backward stepwise selection and the remaining variables were used to cross-validate the final models. Running elastic net (R package glmnet; Friedman et al. 2010) requires tuning of two parameters, α (describes the balance of ridge and lasso) and λ (degree of regularisation; for more details see: Hastie et al. 2009). For each species and bias manipulation, we estimated the best combination of α and λ by 5-fold cross-validation of 11 models (α ranging from zero to one with an increment of 0.1). For each model, the optimal λ value was determined by fitting a series of cross-validated models to a range of λ values (‘regularisation path’; by default 100 values estimated from the data) and the λ value that showed the minimum mean cross-validated error was used for predictions. Similarly, the α value with the lowest error (Poisson deviance) was selected (Supplementary material Appendix 1, Fig. A6). To report the performance of the elastic net (and for comparisons with the results of other techniques), the selected values of α and λ were used to run the cross-validation models manually. Due to the computational limitation of explicitly using a user-defined λ during model training (stated by glmnet help page), models were fitted without providing a λ-value, allowing the fit of many models over the regularisation path. For prediction we used the optimal λ-value estimated from cross-validation. The best-estimated values of α are shown in Supplementary material Appendix 1, Table A2. Lasso (α = 1) rather than ridge was chosen for almost half of the species in environment-only and accessibility models. For the effort model, only three species had ridge models (α = 0) and all others had α-values 1. Maxent default settings were adapted according to advice in the literature (Merow et al. 2013, Radosavljevic and Anderson 2014). Maxent, by default, uses combinations of feature classes (transformations of covariates: linear “L”, quadratic “Q”, hinge “H”, threshold “T”, and product “P”) depending on the number of presence locations available, allowing for complex species-environment relationships (Phillips and Dudík 2008). We adapted functions from the ENMeval package (Muscarella et al. 2014) to run Maxent models on cross-validation at different complexity levels and feature class combinations (48 models = 8 regularisation multiplier values ranging from 0.5 to 4 with increment of 0.5 × 6 feature class combinations [L / LQ / H / LQH / LQHP / LQHPT]). The default number of background points used by Maxent (10,000) is insufficient to represent the environmental variability in large study areas as used here (Renner and Warton 2013), and hence all potential grids were considered as background points. We used clamping while predicting to the left-out fold, meaning that if any value of a covariate is beyond its training range, these values will be replaced by the closest value during training (Anderson and Raza 2010). For each species and bias manipulation, the combination of regularisation multiplier and feature class that shows the highest mean testing-AUC (on cross-validation) were selected for predictions. The optimum combinations of Maxent’s feature classes and regularisation multiplier are shown in Supplementary material Appendix 1, Table A1. The selected feature classes deviated from Maxent’s default: 10 species always had simple features (L/LQ) for all bias manipulations, with the effort model showing more complex features in some cases. Moreover, the best-estimated regularisation multiplier also deviated from the default value of one, with many species having values >1, indicating little overfitting. Model evaluation We evaluated the models using threshold-independent (the area under the ROC curve, AUC) and threshold-dependent (True Skill Statistic, TSS) metrics. We tried to avoid some of the metrics' known caveats through 1) the use of species-specific study areas; 2) block cross-validation minimising environmental extrapolation; and 3) using the same blocking structure to run different modelling techniques and bias manipulations of the same species. In presence-only SDMs, presences comprise only a tiny fraction of the background. During evaluation, the use of too many background points with equal weights for commission and omission error can distort the evaluation (Lobo et al. 2008). To be able to compare AUCs (computed using the dismo R package) among models for the same species, we set the numbers of test presences and random test background (test data prevalence) to a ratio of 1:20, and repeated this 100 times to average stochastic effects. TSS was calculated using the threshold that maximises the sum of sensitivity and specificity (Liu et al. 2013), again using a constant ratio of presences and background. AUC and TSS were calculated for each combination of species, algorithm and bias manipulation. We summarised the overall effect of different bias manipulations in two ways, using linear mixed-effect models (R-package lme4; Bates et al. 2015). Firstly, to evaluate the effect of incorporating bias covariates into the model, we compared evaluation of the environment-only models (without any bias correction) to those of bias-accounted models (accessibility and effort) without conditioning on a particular value of the bias covariate during prediction (modelling evaluation). Second, to explore the effect of bias correction on model evaluation, we performed similar comparisons, with bias-free predictions of the bias-accounted models instead (bias correction evaluation; for details see Supplementary material Appendix 6). We expect models incorporating the impact of bias to outperform (signified by a higher AUC or TSS) those that do not. The problem is that we do not have bias-free data, and therefore our test data exhibit the same bias as the data used to create the model. Under these circumstances, models allowing for bias may be worse than those that do not. In either case, the difference between bias-corrected and control models is a measure of the impact of bias. In the mixed models, we used model evaluation as response variable, species as random effect, and model type, bias correction (and their interactions), total number of training presences, and total number of pixels at per-species study area (range size) as fixed effects (for more details, see Supplementary material Appendix 6). Results Validation by block cross-validation vs. by Egypt hold-out The overall mean AUC is 0.88 ± 0.08 on cross-validation (0.75 ± 0.13 in Egypt’s evaluation; Figs. 2-3). Because results for AUC and TSS were similar, we present only results for AUC (see Supplementary material Appendix 1, Fig. A11 and supplement for TSS). Model evaluations in Egypt were, on average, poorer than on cross-validation. However, both types of evaluation (using ‘bias-corrected’ predictions) show positive correlation (Kendall’s tau; n=21 – Supplementary material Appendix 1, Fig. A7), with highest consistency for the accessibility model. As expected, species with fewer occupied cells tend to have higher evaluation variability (uncertainty), and hence provide less robust analyses (Supplementary material Appendix 1, Figs. A9-10). Effect of bias corrections: mixed-effect model analysis across species and model types Sampling-bias model explained most of the variation in model validation, followed by the modelling algorithm, and their interaction (Fig. 2 and Supplementary material Appendix 6.1). The total number of training presence (positive effect) and the range size (negative effect) were much less important. On average, the accessibility model performed much better than environment-only and effort models (Figs. 2 and Supplementary material Appendix 6.1). For mean cross-validations, effort models also had relatively higher AUCs than the environment-only model (the difference is much smaller for elastic net and GLM; Figs. 2a). However, for validation in Egypt, effort models were similar to the environment-only model (Figs. 2b , see also Supplementary material Appendix 6.1). We can also evaluate the performance of the sampling-bias models when their bias covariates are set to a fixed value, i.e. when predicting to bias-free data. However, since we have no reference survey data to compare them to, these predictions reveal the impact of bias-correction on the model, rather than assess the model’s performance. These evaluations are presented in Fig. 3 and Supplementary material Appendix 6.2. Effect of model type The predictive performance of the three modelling algorithms showed fair to moderately high correlation coefficient (most r > 0.6, all significant; Fig. 4). For mean cross-validation evaluations, GLM performed best, followed by elastic net and Maxent. The order is reversed for independent validations in Egypt: Maxent has highest AUC-values, followed by elastic net, then GLM (Figs. 2-4, and Supplementary material Appendix 6).  Discussion Sampling bias, if not corrected for effectively, can substantially affect the predicted intensity and model evaluation of SDMs. Our results suggest that accessibility bias covariates describe the bias in our focal species’ training data well, compared to sampling efforts, and their use led to higher validation scores. Evaluated on data that are themselves spatially biased, we find that removing sampling bias did not improve the predictive performance on cross-validation or in Egypt (Fig. 3), highlighting the limitation of evaluating bias correction without independent bias-free presence-absence data. Collectively, the three modelling algorithms performed similarly well in cross-validation, though predicting to new sites (Egypt) gave Maxent a slight advantage (Figs. 2b and 3b). Sampling-bias correction using presence-only data Few approaches exist to correct for sampling bias, some not applicable when data are sparse. Spatial filtering (i.e. aggregation to single records within some larger buffer: e.g. Anderson and Raza 2010, Boria et al. 2014) might be wasteful when only few presences exist. Removal of clumped data will reduce training sample size and may remove some of the environmental conditions occupied by the species (depending on the heterogeneity of the landscape, selected distance, and pixel size). Target-group background selection strictly assumes that related species were collected with the same method and equipment (Phillips et al. 2009), bias is similar across species (Warton et al. 2013), and species have the same chance of being recorded in all locations (Yackulic et al. 2013). Moreover, it replaces the observer bias with a (spatial) species-richness bias and can be understood as modelling the likelihood of encountering the target species rather than a non-target species (Warton et al. 2013, Warren et al. 2014). Also, it does not distinguish between areas unsuitable for any of the species and areas of low accessibility (no records, but potentially suitable: Warren et al. 2014, Fithian et al. 2015). The target-group background cannot be used when data on the related species are similarly limited, as this increases the risk of extrapolating in environmental space (Mateo et al. 2010, Merow et al. 2013). Bias layers attempt to describe with which biases presence locations were recorded (Phillips et al. 2009, Warton et al. 2013). Our external sampling-effort model (Fig. 1 bottom right) is such an attempt to derive a bias layer (e.g. Elith et al. 2010). The bias layer is currently only implemented in Maxent, and we are not aware of any study that applied a similar approach using other modeling algorithms. To maintain consistency of the analyses across different modelling algorithms we did not consider using a bias layer here. Both bias layer and target-group are used to sample background with the same bias as presence locations, which is a sensitive tuning step that strongly influences model evaluation (Chefaoui and Lobo 2008). We think that model-based bias correction, as applied here, is more plausible in data deficient situations and if applied efficiently, it frees us from artificially manipulating presences or background points, and rather focus on sampling bias correction. The effectiveness of bias correction depends on whether bias covariates are actually able to describe the bias in the available data. Our use of either accessibility or sampling effort bias covariates did not lead to notably different conclusions (see Fig. 5). Using sampling-effort bias covariate led to, on average, little to moderate evaluation changes (Figs. 2-3). It is in fact a model-based version of the target-group background approach, without strict selection of the background data points. The assumption that closely related species have similar bias may not hold. Surprisingly, most of the available presences of our focal species were located in areas of low to moderate estimated efforts (Supplementary material Appendix 1, Fig. A13), but still strongly biased towards roads and cities (and, to less extent, protected areas; Supplementary material Appendix 1, Fig. A13). Ideally, evaluating bias correction requires independent ‘bias-free’ presence-absence testing data (Phillips et al. 2009). Such data are typically not available, especially in developing countries, and removing bias from test data is very difficult (Dudík et al. 2005, Smith 2013). On average, validation of bias-corrected predictions using well-structured, independent presence-absence data from rigorous surveys led to improved predictions (Elith and Leathwick 2007, Phillips et al. 2009, Mateo et al. 2010, Syfert et al. 2013, Warton et al. 2013, Boria et al. 2014). However, this improvement is not happening in all situations and depends on the modelling conditions, species prevalence, validity of assumptions, and how effective bias covariates are in describing bias in training presences (Phillips et al. 2009, Warton et al. 2013, Fourcade et al. 2014). As covariates that affect sampling may also affect the distribution of a species (e.g. avoiding deserts), no method can fully correct for sampling bias in presence-only data without affecting the niche model (Merow et al. 2014, Guillera-Arroita et al. 2015). Correction of sampling bias leads to larger areas of suitable habitats due to higher suitability estimates in low-accessible sites (e.g. Phillips et al. 2009, Warton et al. 2013; Fig. 5). Predictions at such sites are of lower reliability and should be interpreted with caution (Supplementary material Appendix 1, Fig. A16 for an example). They can be used to guide future surveys and conservation planning, but not for taking serious conservation decisions (Guisan et al. 2006). Evaluations using spatial-block cross-validation vs. Egyptian hold-out We used spatial-block cross-validation to avoid overestimating model performance and underestimating predictive errors (Bahn and McGill 2013, Renner et al. 2015, Roberts et al. 2017). Block cross-valuations suggested a better model performance than evaluation on Egyptian data. This can be explained by differences in sample size and by environmental variability: on average, our cross-validation models had a larger mean number of testing presences and higher environmental variability compared to evaluations at smaller scale (Table 1). However, evaluations at both scales were positively correlated, which supports the idea that evaluations on cross-validation (larger extent) can be indicative of performance at the local scale. Evaluation metrics & predictive performance Using AUC and TSS led to consistent conclusions. The overall evaluation scores are fair to high, taking into account that spatially independent evaluation yields lower scores compared to commonly used random split (Radosavljevic and Anderson 2014). In presence-only SDMs, the use of a particular value for defining good models become unreliable (Yackulic et al. 2013), due to higher uncertainty of estimates at background points and as the number of background points used is not fixed. In this study, we maintained constant test-data prevalence (1:20) across all comparisons, which gave very similar results, across all species, to the standard default approach of computing AUC (see Supplementary material Appendix 1, Fig. A14). For some species, the spread over the 100 repetitions was very noticeable, however, making our approach somewhat more robust. The improved predictive performance reported by Syfert et al. (2013), for example, may be in part attributable to much higher testing prevalence. Also in our study, the larger the number of available presences and the smaller the study area was, the higher were the predictions scores. Comparison of modelling algorithms The three modelling algorithms applied did not lead to different conclusions and their evaluations were highly correlated. Elastic-net models showed intermediate performance in all situations. GLM (more specifically: the down-weighted Poisson regression (DWPR) with variable selection) had the highest evaluation on cross-validation, which may suggest that GLM-DWPR is more powerful with a larger number of testing presences. However, GLMs predictions in Egypt were less nuanced compared to other two modelling algorithms (Supplementary material Appendix 1, Fig. A15a, for example), which explains why it ranked lowest on the Egyptian hold-out. Maxent had the lowest prediction error on the Egyptian data, which suggests its transferability to situations of low extrapolation. However, Maxent had the lowest evaluation on cross-validation, possibly due to clamping. We applied ‘clamping’ to constrain the response beyond the training range, which changes some of the predicted values. Clamping can affect model evaluation (the ranking of the predicted values) and, while its effect has not been well-explored in the literature, it seems to depend on the shape of the response curve (at both ends), the importance of the covariate, and how much environmental extrapolation occurs. We expect the effect of clamping to be small in the Egyptian hold-out compared to the cross-validation, due to less extrapolation. When correcting for bias, no interaction should be allowed between bias and other covariates so that, for prediction, the bias can be corrected for without affecting the other covariates (Warton et al. 2013). For GLM and elastic net, we have full control of the models’ interactions. However, the version of Maxent used here (3.3.3k) does not enable users to select which interactions to use, and the use of the product feature inevitably enables all pairwise interactions. In further applications, it may be recommended to disable the ‘product’ feature class while correcting for sampling bias. However, its use here did not lead to different conclusions compared to GLM and elastic net, suggesting that bias-suitability interactions were of limited importance. During the reviewing of this manuscript, an open-source version of Maxent (maxnet R package: Phillips et al. 2017) was released that uses the glmnet package for L1-regularization. Maxnet provides flexibility for specifying the interactions to be used, so it is possible to exclude interactions between environmental and bias variables (making them similar to our elastic-net models, but with more flexibility for other feature classes implemented in Maxent, e.g. the hinge). This early version of maxnet implements an Infinitely-Weighted Logistic Regression (IWLR; Fithian and Hastie 2013) and only L1-regularization (lasso); however, further extensions are possible, e.g. the implementation of DWPR and elastic net (similar to our elastic-net models). We implemented the down-weighted Poisson regression using both GLM and elastic net. As Poisson models, their predictions have no upper bound and may thus yield extreme predictions. We reported the existence of a few extreme predicted intensities for many species, which makes it difficult to plot their predictions on a linear scale (see Fig. 5 and Supplementary material Appendix 1, Fig. A15b for an example). In contrast, Maxent puts a constraint on the moments of the predictions, making them less subject to extreme values (Phillips et al. 2006). Conclusion Data-sparse regions pose challenges to modelling species distributions, exacerbated by noticeable sampling biases. We recommend the use of model-based bias correction in data-sparse situations, in which other bias corrections methods are not possible; however, the effectiveness of bias correction depends on whether bias covariates are actually able to describe the bias in the available data. Using covariates to describe site accessibility improved prediction to spatially independent hold-outs, compared to environment-only or effort models. Augmenting local records with data from across the species’ range allowed us to make consistently high-quality predictions to hold-out data from an entire country (in this case Egypt). Bias-free predictions can enhance future conservation planning and target future surveys when limited resources are available to cover large study areas. However, due to possible lower certainty at unsurveyed locations, they should be used cautiously (maps including bias are of use only during model cross-validation). Without survey-based presence-absence data, no complete evaluation of the quality of bias corrections can be attempted. Down-weighted Poisson regression as well as the statistically equivalent Maxent approach led to similar results, with more flexibility in the elastic-net models (e.g. degree of shrinkage, question-led specification of non-linear effects and interactions). More important than the specific algorithm is to use the point-process modelling framework as such. Acknowledgements We would like to express our sincere thanks for Petr Benda (Charles University) for comments on the distribution of the bat species. An earlier version of the manuscript was improved by comments of Prof. Francis Gilbert (University of Nottingham). We are also grateful to Belarmain Fandohan, David R. Roberts, and Simone Ciuti for long discussions and comments on the manuscript. This manuscript benefitted from comments of the subject editor Robert P. Anderson and two anonymous reviewers. A. El-G. is sponsored by the German Academic Exchange Service (DAAD) through a GERLS scholarship. This work was partially performed on the computational resource ‘bwUniCluster’ funded by the Ministry of Science, Research and Arts and the Universities of the State of Baden-Württemberg, Germany, within the framework program bwHPC. References Anderson, R. P. 2003. Real vs. artefactual absences in species distributions: Tests for Oryzomys albigularis (Rodentia: Muridae) in Venezuela. – Journal of Biogeography 30: 591–605. Anderson, R. P. and Raza, A. 2010. The effect of the extent of the study region on GIS models of species geographic distributions and estimates of niche evolution: preliminary tests with montane rodents (genus Nephelomys) in Venezuela. – J Biogeogr 37: 1378–1393. Bahn, V. and McGill, B. J. 2013. Testing the predictive performance of distribution models. – Oikos 122: 321–331. Barve, N. et al. 2011. The crucial role of the accessible area in ecological niche modeling and species distribution modeling. – Ecol Model 222: 1810–1819. Basuony, M. I. et al. 2010. Mammals of Egypt: atlas, red data listing & conservation. BioMAP & CultNat, EEAA & Bibliotheca Alexandrina, Cairo. 286 pp. Bates, D. et al. 2015. Fitting Linear Mixed-Effects Models Using lme4. – J Stat Softw 67. Boria, R. A. et al. 2014. Spatial filtering to reduce sampling bias can improve the performance of ecological niche models. – Ecol Model 275: 73–77. Chefaoui, R. M. and Lobo, J. M. 2008. Assessing the effects of pseudo-absences on predictive distribution model performance. – Ecol Model 210: 478–486. Dudík, M. et al. 2005. Correcting sample selection bias in maximum entropy density estimation. Advances in Neural Information Processing Systems 18. The MIT Press, pp. 323-330. URL: http://papers.nips.cc/paper/2929-correcting-sample-selection-bias-in-maximum-entropy-density-estimation.pdf. Elith, J. et al. 2010. The art of modelling range-shifting species. – Methods Ecol Evol 1: 330–342. Elith, J. et al. 2011. A statistical explanation of MaxEnt for ecologists. – Divers Distrib 17: 43–57. Elith, J. and Leathwick, J. 2007. Predicting species distributions from museum and herbarium records using multiresponse models fitted with multivariate adaptive regression splines. – Divers Distrib 13: 265–275. Fithian, W. et al. 2015. Bias correction in species distribution models: pooling survey and collection data for multiple species. – Methods Ecol Evol 6: 424–438. Fithian, W. and Hastie, T. 2013. Finite-sample equivalence in statistical models for presence-only data. – Ann Appl Stat 7: 1917–1939. Fourcade, Y. et al. 2014. Mapping species distributions with MAXENT using a geographically biased sample of presence data: a performance assessment of methods for correcting sampling bias. – PloS one 9: e97122. Friedman, J. H. et al. 2010. Regularization paths for generalized linear models via coordinate descent. – Journal of statistical Software 33 - http://www.jstatsoft.org/v33/i01/paper. Gaiji, S. et al. 2013. Content assessment of the primary biodiversity data published through GBIF network: status, challenges and potentials. – Biodiversity Informatics 8: 94–122. Guillera-Arroita, G. et al. 2015. Is my species distribution model fit for purpose? Matching data and models to applications. – Global Ecol Biogeogr 24: 276–292. Guisan, A. et al. 2006. Using Niche-Based Models to Improve the Sampling of Rare Species. – Conservation Biology 20: 501–511. Hastie, T. et al. 2009. The elements of statistical learning: data mining, inference, and prediction - Springer New York. Liu, C. et al. 2013. Selecting thresholds for the prediction of species occurrence with presence-only data. – J Biogeogr 40: 778–789. Lobo, J. M. et al. 2008. AUC: a misleading measure of the performance of predictive distribution models. – Global Ecol Biogeogr 17: 145–151. Mateo, R. G. et al. 2010. Profile or group discriminative techniques? Generating reliable species distribution models using pseudo-absences and target-group absences from natural history collections. – Divers Distrib 16: 84–94. Merow, C. et al. 2013. A practical guide to MaxEnt for modeling species’ distributions: what it does, and why inputs and settings matter. – Ecography 36: 1058–1069. Merow, C. et al. 2014. What do we gain from simplicity versus complexity in species distribution models? – Ecography 37: 1267–1281. Muscarella, R. et al. 2014. ENMeval: an R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. – Methods Ecol Evol 5: 1198–1205. Newbold, T. 2010. Applications and limitations of museum data for conservation and ecology, with particular attention to species distribution models. – Progress in Physical Geography 34: 3–22. Pearce, J. L. and Boyce, M. S. 2006. Modelling distribution and abundance with presence-only data. – J Appl Ecol 43: 405–412. Phillips, S. J. et al. 2006. Maximum entropy modeling of species geographic distributions. – Ecol Model 190: 231–259. Phillips, S. J. et al. 2009. Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data. – Ecological Applications 19: 181–197. Phillips, S. J. et al. 2017. Opening the black box: an open-source release of Maxent. – Ecography. Phillips, S. J. and Dudík, M. 2008. Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. – Ecography 31: 161–175. Ponder, W. F. et al. 2001. Evaluation of Museum Collection Data for Use in Biodiversity Assessment. – Conservation Biology 15: 648–657. Radosavljevic, A. and Anderson, R. P. 2014. Making better Maxent models of species distributions: complexity, overfitting and evaluation. – J. Biogeogr. 41: 629–643. Renner, I. W. et al. 2015. Point process models for presence-only analysis. – Methods Ecol Evol 6: 366–379. Renner, I. W. and Warton, D. I. 2013. Equivalence of MAXENT and Poisson point process models for species distribution modeling in ecology. – Biometrics 69: 274–281. Roberts, D. R. et al. 2017. Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure. – Ecography 28. Smith, A. B. 2013. On evaluating species distribution models with random background sites in place of absences when test presences disproportionately sample suitable habitat. – Diversity Distrib. 19: 867–872. Stolar, J. and Nielsen, S. E. 2015. Accounting for spatially biased sampling effort in presence-only species distribution modelling. – Divers Distrib 21: 595–608. Syfert, M. M. et al. 2013. The effects of sampling bias and model complexity on the predictive performance of MaxEnt species distribution models. – PloS one 8: e55158. Warren, D. L. et al. 2014. Incorporating model complexity and spatial sampling bias into ecological niche models of climate change risks faced by 90 California vertebrate species of concern. – Divers Distrib 20: 334–343. Warton, D. and Aarts, G. 2013. Advancing our thinking in presence-only and used-available analysis. – The Journal of Animal Ecology 82: 1125–1134. Warton, D. I. et al. 2013. Model-based control of observer bias for the analysis of presence-only data in ecology. – PloS one 8: e79168. Yackulic, C. B. et al. 2013. Presence-only modelling using MAXENT: when can we trust the inferences? – Methods Ecol Evol 4: 236–243. Table 1: List of Egyptian bat species used in this study, with total number of records available, number of unique records (in parentheses: the number of records used to train the model on cross-validation/those kept aside for independent evaluation, i.e. ‘Egyptian records’), number of occupied grid cells at the resolution of 5 × 5 km (number of cells outside/inside Egypt), and best estimated number of background points used to run the DWPR-GLM and elastic-net models. See main text and Supplementary material Appendix 1, Fig. A5 for more details. Species # Records (Total) # Records (unique) # occupied cells # background points (×1000) Asellia tridens (Trident Leaf-nosed Bat) Barbastella leucomelas (Sinai Barbastelle) Eptesicus bottae (Botta’s Serotine Bat) Hypsugo ariel (Fairy Pipistrelle) Nycteris thebaica (Egyptian Slit-faced Bat) Nycticeinops schlieffeni (Schlieffen's Bat) Otonycteris hemprichii (Hemprich’s Long-eared Bat) Pipistrellus deserti (Desert Pipistrelle) Pipistrellus kuhlii (Kuhl's Pipistrelle) Pipistrellus rueppellii (Rueppell’s Pipistrelle) Plecotus christii (Desert Long-eared Bat) Rhinolophus clivosus (Arabian Horseshoe Bat) Rhinolophus hipposideros (Lesser Horseshoe Bat) Rhinolophus mehelyi (Mehely’s Horseshoe Bat) Rhinopoma cystops (Lesser Mouse-tailed Bat) Rhinopoma microphyllum (Greater Mouse-tailed Bat) Rousettus aegyptiacus (Egyptian Fruit Bat) Tadarida aegyptiaca (Egyptian Free-tailed Bat) Tadarida teniotis (European Free-tailed Bat) Taphozous nudiventris (Naked-Bellied Tomb Bat) Taphozous perforatus (Tomb Bat) Figure 1: Flowchart of analyses in this study, illustrated with data for Asellia tridens. Sampling-bias models and modelling algorithms were combined factorially. Only results for validation with AUC are presented in the manuscript, while TSS-results are given in the appendix. Figure 2: Mean AUC of each species calculated either on 5-fold spatial block cross-validation (a) or in Egypt (b). Each plot compares evaluation of environment-only models to bias-accounting models (effort and accessibility), without correcting for sampling bias (Modelling evaluation; see Appendix 5.1). Each species is represented by different symbols (see Supplementary material Appendix 1, Fig. A7 for species names). Red lines indicate overall mean AUC for each modelling algorithm and bias manipulation applied. For evaluations of bias-free prediction, see Fig. 3. Figure 3: Species mean AUC calculated either on cross-validation (a) or in Egypt (b), after sampling bias correction (using bias-free prediction; for details, see Supplementary material Appendix 6.2). Each species is represented by different symbols (see Supplementary material Appendix 1, Fig. A7 for species names). Red lines indicate the overall mean AUC at each modelling algorithm and bias models applied. Figure 4: Kendall’s correlation of the per-species mean AUC between pairs of modelling algorithms. Each species is represented by different symbols (see Supplementary material Appendix 1, Fig. A7 for species names) with colours referring to bias models (using predictions of environment-only model and bias-free prediction of accessibility and effort model). ‘M’ indicates the overall mean evaluation. Top row panels are mean evaluation using spatial block cross-validation, while those in the bottom row are independent evaluations in Egypt. Figure 5: Mean cross-validated predicted distribution of Otonycteris hemprichii, of different modelling algorithms (rows) and bias models (columns). Maps were rescaled to relative intensities between zero and one, as different modelling algorithms do not have the same scale. Darker colour indicates higher predicted relative intensity. Blue points (in the top left panel) represents available records used for cross-evaluation model training. For visualisation, extreme values (> 0.9995 quantile of predicted values) were replaced with their next smaller value, as GLM and elastic net are subject to some few extreme predictions (see main text, and Supplementary material Appendix 1, Fig A15B for a comparison with unlimited predictions). For predicted maps in Egypt see Supplementary material Appendix 1, Fig. A15A.

Wrong, but useful: regional species distribution models may not be improved by range-wide data under biased sampling
Ahmed El-Gabbas
Carsten F. Dormann
https://doi.org/10.1002/ece3.3834
10.1002/ece3.3834
Species distribution modeling (SDM) is an essential method in ecology and conservation. SDMs are often calibrated within one country's borders, typically along a limited environmental gradient with biased and incomplete data, making the quality of these models questionable. In this study, we evaluated how adequate are national presence-only data for calibrating regional SDMs. We trained SDMs for Egyptian bat species at two different scales: only within Egypt and at a species-specific global extent. We used two modeling algorithms: Maxent and elastic net, both under the point-process modeling framework. For each modeling algorithm, we measured the congruence of the predictions of global and regional models for Egypt, assuming that the lower the congruence, the lower the appropriateness of the Egyptian dataset to describe the species' niche. We inspected the effect of incorporating predictions from global models as additional predictor (“prior”) to regional models, and quantified the improvement in terms of AUC and the congruence between regional models run with and without priors. Moreover, we analyzed predictive performance improvements after correction for sampling bias at both scales. On average, predictions from global and regional models in Egypt only weakly concur. Collectively, the use of priors did not lead to much improvement: similar AUC and high congruence between regional models calibrated with and without priors. Correction for sampling bias led to higher model performance, whatever prior used, making the use of priors less pronounced. Under biased and incomplete sampling, the use of global bats data did not improve regional model performance. Without enough bias-free regional data, we cannot objectively identify the actual improvement of regional models after incorporating information from the global niche. However, we still believe in great potential for global model predictions to guide future surveys and improve regional sampling in data-poor regions.
Species distribution models (SDMs) are statistical methods that relate species information (either presence-only or presence–absence) to environmental variables to infer spatially explicit habitat suitability. They are being used intensively as a standard tool for estimating potential range shifts under climate change, assessing invasion risk, locate future survey sites, and conservation planning and prioritization (Araújo, Alagador, Cabeza, Nogués-Bravo, & Thuiller, 2011; Guisan & Zimmermann, 2000; Guisan et al., 2013; Rodríguez, Brotons, Bustamante, & Seoane, 2007; Thuiller et al., 2005). Although these methods have limitations and uncertainties (Araújo & Guisan, 2006; Dormann, Purschke, Márquez, Lautenbach, & Schröder, 2008; Guisan & Thuiller, 2005), they constitute the best available tools when not much detailed information on the ecology and physiology of the species is available (Warren, Wright, Seifert, Shaffer, & Franklin, 2014).
In developing countries, the majority of species sightings are scattered, opportunistic, and recorded mainly in museum catalogues, personal collections, and the literature. Due to political instability and limited funds dedicated to wildlife conservation (amongst other reasons), there is no systematic nation-wide sampling scheme for collecting biological information in most developing countries. Many of these countries do not share their biodiversity data, making them highly under-represented at international data depositories, such as the Global Biodiversity Information Facility (GBIF), with many more records from countries with high GDP (Newbold, 2010). Furthermore, data from developing countries are particularly (but not exclusively) spatially biased (more records from accessible locations near roads and cities) and taxonomically biased (toward larger or charismatic species). Spatial bias poses a problem for SDMs, which, in their default approach, assume that available presence locations represent a random (representative) sample in the environmental/geographical space, with no spatial dependencies (Elith et al., 2011; Renner et al., 2015). This assumption is hardly ever met due to sampling bias, imperfect detectability and spatial auto-correlation (Guillera-Arroita et al., 2015). When high sampling bias exists, SDM predictions provide an estimate not necessarily of the species suitability, but more of the patterns of the sampling effort and detectability (Elith et al., 2011; Yackulic et al., 2013). Several methods have been proposed to correct for sampling bias (e.g., target-group background: Phillips et al., 2009; spatial filtering: Anderson & Raza, 2010; sampling bias predictors: Warton, Renner, & Ramp, 2013); however, no method seems to be able to fully correct for sampling bias in presence-only data (El-Gabbas & Dormann, 2017; Merow et al., 2014).
One of the major challenges of SDM studies is how to determine the extent of the study area appropriately. Study area should be objectively determined to cover accessible areas by the species within its known complete range, allowing for wider range of environmental variation and extremes occupied by the species (Barve et al., 2011; Raes, 2012; Sánchez-Fernández, Lobo, & Hernández-Manrique, 2011). However, it is common that study areas are unjustifiably determined based on geographical or political borders for regional/local conservation actions, resulting in models calibrated with a limited range of environmental conditions that do not capture much of the species' niche and hence is insufficient to describe its environmental tolerance (Raes, 2012; Titeux et al., 2017). This leads to the truncation of the estimated response curves, underrepresentation of areas of suitable habitats, and limiting the predictive power of the models (Sánchez-Fernández et al., 2011; Thuiller, Brotons, Araújo, & Lavorel, 2004). This is more problematic when the aim of the study is to extrapolate beyond the training range, either in time or space (Barbet-Massin, Thuiller, & Jiguet, 2010; Thuiller et al., 2004), or in situations where available data are few, opportunistic, or with high (typically unknown) sampling bias. The paucity of available records in developing countries, coupled with clear signs of sampling bias and limited local environmental gradients, makes it challenging to establish robust SDMs for a variety of taxonomic groups at the national scale.
In this article, we evaluate the adequacy of regional presence-only data (in this case from within a developing country's political borders) for constructing SDMs. More specifically, we compare bat occurrence predictions from regional and global SDMs for the country of Egypt, in many respects exemplary for developing countries. Egypt shows much lower environmental variability compared to the global extents of the species (see Figures 1 and S1) and comprises only a small proportion of available global records. This makes the quality of regional SDMs, that is, those built only on the sparse Egyptian data, questionable. Global models (at species-specific global range) should in this case be more reliable than regional models (in Egypt) in describing the climatic niche of species because they are calibrated with a much higher number of presences and capture a much wider range of occupied (or, more generally, accessible) environmental conditions (Pearson, Dawson, & Liu, 2004). Thus, we evaluate predictions from regional and global SDMs for Egypt, arguing that the less similar they are, the more the local data describe sampling effort rather than the ecology of bats. Furthermore, we investigate how much correction for sampling bias (using bias predictors, in both regional and global SDMs) helps to improve the local predictions for Egypt.
Predictions from global models interpolated to Egypt represent a spatial-explicit information on the species potential distribution that is independent from regional data available from Egypt, and thus can be useful to improve predictions of regional models when used as additional predictors (cf. “informative offset”: Merow, Allen, Aiello-Lammens, & Silander, 2016). We explore how much global predictions (interpolated to Egypt) improve Egyptian regional models when used as predictor “prior” to describe the environmental niche (again, with and without correcting for sampling bias).
Study design and species
This study builds on a comparison of methods to correct for sampling biases (El-Gabbas & Dormann, 2017), adding an evaluation of regional species distribution models based on national records. We collected records for Egyptian bat species (from within Egypt and their global extents) from different sources (Appendix S1 and El-Gabbas & Dormann, 2017). Four species with fewer than eight unique sightings in Egypt were excluded from the analyses, yielding a total of 17 species (Table S1). For the selected species, we created regional models using presence locations and environmental data only for Egypt (“regional SDMs”). “Regional” refers here to a geographic extent much smaller than the range of the species, but of much coarser grain than a local dataset. We also created analogous models across the global range (“global SDMs”): These models were made for each species-specific global extent (a buffered bounding box around all global records), excluding Egyptian records to maintain independence (and to allow for valid comparisons) between the predictions of the regional and global models (see below; and El-Gabbas & Dormann, 2017 for details). Both scales are nested in geographical and environmental space: Our regional models are calibrated within a subset of each species-specific global extent. At either scale, we used two modeling algorithms under the point-process modeling framework (Maxent and elastic net; Renner et al., 2015), with two options on dealing with sampling bias (with and without bias correction), and evaluated the results using spatial-block cross-validation (Roberts et al., 2017).
Environmental variables
Potential environmental predictors (at the total study area covering both scales) and species records were projected into Mollweide equal-area projection at a resolution of 5 × 5 km2. Using the same pixel size and projection maintains consistency of the analyses between regional and global models (Budic, Didenko, & Dormann, 2016). As the correlation between predictors varies from one study area to another, different environmental predictor combinations were used at regional and global scales. Some predictors were not useful at the regional scale, and hence were excluded a priori; for example, precipitation of driest month does not show any variability across Egypt because most of Egypt receives no precipitation at all in summer, reflecting its hyper-arid climate (El-Gabbas, Baha El Din, Zalat, & Gilbert, 2016). We ensured minimum multi-collinearity at both scales by selecting only predictors that maintain a maximum generalized variance inflation factor value less than 3 (see Table S2 for the list of predictors used at either scale).
Modeling algorithms
We used two modeling algorithms: Maxent and elastic net. Maxent (Phillips & Dudík, 2008; v3.3.3k) is a machine-learning presence-background SDM algorithm. It outperforms other presence-only SDM algorithms, especially at smaller sample sizes (e.g., Wisz et al., 2008), due to its use of (some form of) lasso regularization. Elastic net (Friedman, Hastie, & Tibshiani, 2010) is an extension of GLMs that uses “lasso” and “ridge” regularization rather than AIC to select the most suitable model, and hence is similarly resistant to overfitting. We applied both algorithms under the point-process modeling framework following recommendations of Renner et al. (2015), changing some of Maxent's default settings (e.g., to “noautofeature,” “noaddsamplestobackground,” and “noremoveduplicates”), and used the implementation of “down-weighted Poisson regression” for elastic-net models. For each calibrated model of either algorithm, we adjusted against unnecessary complexity (Merow et al., 2014) using five-fold spatial-block cross-validation, estimating the best combination of Maxent's feature classes and regularization multiplier based on maximizing the mean testing AUC (Muscarella et al., 2014), and the optimum α (which describes the balance between ridge and lasso) for elastic net.
Adjusting for sampling bias
In addition to “environment-only” models (without bias correction), we use two different methods of predicting from models that incorporate bias: “bias-predictor” and “bias-corrected.” In both methods, we use sampling bias predictors as our estimate of bias: three layers describing distances to main roads, cities, and protected areas (Warton et al., 2013). Bias-predictor models use the bias layers simply as an extra set of predictors, and during prediction also their values change. Bias-corrected models try to factor out the bias by setting the bias variables to zero (see Warton et al., 2013). The three options for sampling bias (none, predictor, and correction) were applied to regional and global models, with bias predictors nested for regional scale within the global scale.
Model evaluation and the use of spatial priors
We evaluated regional model performance using AUC as a threshold-independent metric. Despite the criticism of the use of AUC to evaluate the performance of presence-only SDMs (e.g., Lobo, Jiménez-Valverde, & Real, 2008), our use of AUC for comparisons between models of the same species, predictors, and study area is valid (Anderson & Gonzalez, 2011; Wisz et al., 2008). We did not use AUC to quantify model performance (goodness of fit), but rather as a measure of the relative ranking of predictions at testing presence and background locations. We calculated AUC on five-fold spatial-block cross-validation to maintain spatial independence between training and testing data (Fithian, Elith, Hastie, Keith, & O'Hara, 2015; Roberts et al., 2017): The same blocking structure (how spatial blocks are distributed into cross-validation folds) is used for each species, with balanced prevalence among blocks and same block sizes, allowing for valid AUC comparisons for the same species. The mean value of testing AUC on cross-validation is reported.
To quantify the efficacy of Egyptian data to construct SDMs, we calculated the geographical congruence (Schoener's D; Schoener, 1968; Warren, Glor, & Turelli, 2010) between continuous predictions of the global and regional SDMs for Egypt (scaled to sum to one; without and with bias correction). Our assumption is that the higher the geographical congruence, the more suitable the Egyptian records are to parameterize regional models. When assessing the congruence between maps we used all three bias options, while for regional comparisons based on AUC we only used the first two models (environment-only and bias-predictor), due to the lack of bias-free testing-data from Egypt required to evaluate bias-corrected predictions. Geographical congruence and AUC gave similar results, indicating that geographical congruence indeed measured how similarly well, not how similarly poorly models predicted.
We then measured the improvement of regional SDMs after incorporating a spatial-explicit information on the global climatic niche. More specifically, for each species we used predictions from the global SDM interpolated to Egypt (i.e., not using the Egyptian data, and thus referred to hereafter as “prior”) as an additional predictor to create a new set of regional models. We had three types of priors representing the predictions of global models for Egypt: 1) from the environment-only model, “Priorenv-only”; 2) a prediction incorporating the bias layer as a predictor to adjust for sampling bias, “Priorbias-predicted”; and 3) a prediction from a model that has factored out bias, “Priorbias-corrected”. Modeling algorithms were not mixed, that is, global models from Maxent were used only for regional models with Maxent, and analogously for elastic-net models. We quantified the improvement due to priors in two ways. First, we measured changes in model performance (AUC). Secondly, we calculated the map congruence between regional models' predictions in Egypt with and without incorporating priors: the higher the map congruence, the lower the contribution of the prior to the regional SDM. One-tailed paired t-test (df = 16) was used for comparisons between each pair of modeling algorithms, sampling bias options, and changes in AUC and map congruence.
The relative importance of environmental variables (permutation importance calculated by Maxent) varied at global and regional scales. When incorporated, the accessibility bias predictors at both scales had high Maxent permutation importance (particularly, “distance to cities” was of significantly higher importance than all but one variable [p .05; nonsignificant only for Bio4 at global scale and Bio6 at regional scale], and “distance to roads” which had a significantly higher average importance than three different environmental variables at either scales; Figure 2). Furthermore, the response of species to environmental predictors was, unsurprisingly, different at both scales. For example, for Eptesicus bottae at the global scale, the response to precipitation of the coldest quarter increased sharply at low precipitation values (approx. 0–130 mm), then remained high or decayed depending on whether the global bias predictors were used or not, respectively (Figure S2a). At the regional scale, however, the species response was highest at extremely low precipitation values (around 10 mm), then declined sharply (Figure S2c).
Global versus regional SDMs
Different areas were identified as suitable in models either using data from the full range or just from Egypt, with low geographic congruence between the predictions of global and regional models for Egypt (Figure 3). The incorporation of bias predictors (at both scales) did not lead to substantial congruence improvement (yet statistically significant; all p .01). The congruence was highest when bias-corrected models were used (statistically higher than environment-only and bias-predicted models for Maxent and elastic net, p .001). Maxent and elastic net yielded similar values for congruence, with an advantage of Maxent for bias-predictor models (p .05).
The use of prior information from the entire range
The use of priors did not lead to AUC improvement, except when using Priorbias-predictor (p .05; Figure 4a). Results were similar for both Maxent and elastic net, with higher AUC values for Maxent (all p .01). Maxent showed relatively low permutation importance of the different prior variables, except for Priorbias-predictor which had high contributions to the models (all p .0001, although also with high variability; Figure S3, left panel).
The incorporation of prior variables as predictors yielded high geographical congruence between the predictions of regional models without and with priors (Figure 5). However, the congruence values depended on the prior used. The use of Priorenv-only or Priorbias-corrected led to high congruence, indicating little additional information provided by the priors. In contrast, when Priorbias-predictor was used, geographical congruence was less pronounced (p .001), suggesting that here information different from the regional data entered the model. Both Maxent and elastic net produced similar values for congruence, with slightly higher values for elastic net when Priorbias-predictor was used (marginally significant; p = .042).
Correction of regional sampling bias
When regional bias predictors were incorporated into the SDMs, the regional models performed better (higher AUC; all p .05), leading to a negligible effect of priors (Figure 4b). Maxent has relatively higher AUC scores than elastic net (all p .01). However, Priorbias-predictor showed equivalently high AUC values whether or not regional bias predictors were included (p > .7; see Figure 4a,b for a comparison). This was also evident by the much lower permutation importance of prior predictors when regional bias predictors were incorporated, with relatively higher importance for Priorbias-predictor (all p .05; Figure S3, right panel).
Incorporating regional bias predictors led to similar patterns of congruence (between predictions of regional SDMs created with or without priors) to those which did not incorporate bias (Figure 5 vs. Figure S4, light gray boxes), with relatively lower congruence when Priorbias-predictor was used. However, bias-correction (factoring out the bias) did not affect congruence for Maxent, while much lower congruence values were observed for elastic net whichever priors were used (Figure S4, dark gray boxes). In other words, regional bias correction led to less agreement between regional model predictions (with and without priors) for elastic net, regardless of which prior variables were used.
In this study, we evaluated how much improvement to the regional SDMs for Egypt occurs by incorporating additional information (the “priors”) representing the global climatic niche from outside Egypt. First, without providing information on regional bias (no regional bias correction), Priorenv-only and Priorbias-corrected did not lead to improvements in the regional models: Similar AUC values (Figure 4a) and high geographical congruence (Figure 5) imply that they do not provide new information to the regional models. However, the use of Priorbias-predictor led on average to higher AUC and lower geographical congruence, signaling that new information was provided to the models. This was supported in Maxent models by the higher permutation importance of Priorbias-predictor, compared to the other two options of priors (Figure S3, left panel). On the other hand, when regional bias predictors were incorporated, all models had improved AUC, whether or not priors were used (Figure 4b). Regional bias predictors describe the local bias existing in the Egyptian dataset, and their use led to higher AUC, in accordance with other studies (El-Gabbas & Dormann, 2017; Warton et al., 2013). The use of regional bias predictors makes the contribution of priors negligible: Priorenv-only and Priorbias-corrected had an extremely low contribution to these models, only slightly higher for Priorbias-predictor (Figure S3, right panel). Generally, Maxent and elastic net led to very similar results, with slightly higher discrimination ability for Maxent.
Priorbias-predictor implicitly contains information on the regional bias of the records in Egypt, because it represents predictions of equivalent global models calibrated with accessibility bias variables (regional bias variables represent a narrower range than their equivalent variables at global scale). In contrast to bias-free predictions, the use of bias variables as predictors gives higher predicted suitabilities at locations of high accessibility (e.g., closer to roads and cities), which is the reason for high AUC scores when evaluation datasets are similarly biased (Warton et al., 2013). The available dataset for Egyptian bats is spatially-biased, with more records collected near roads and cities (El-Gabbas & Dormann, 2017), and hence Priorbias-predictor describes the available data better than the other two priors. The relatively modest contribution of Priorbias-predictor, and even lower contribution of the other two priors, can be understood as the result of the unavailability of complete, bias-free data from Egypt (see below). Furthermore, Priorenv-only and Priorbias-corrected are highly correlated with some other environmental variables in Egypt (higher than for Priorbias-predictor), particularly for Bio19 (precipitation of coldest quarter) and Bio9 (mean temperature of driest quarter; Figure S5), and hence to a large extent provide redundant information.
The three prior suitabilities show low geographical congruence with their corresponding regional predictions in Egypt (Figures 3 and S6, e.g., maps), meaning they (global models) identify different sites as suitable than do models based on Egyptian records. This can be explained by factors related to model misspecification (e.g., the variables used and violation of model assumptions), the difficulty of modeling widespread species with high accuracy (Stockwell & Peterson, 2002), the low quality of available data, or species-specific reasons (e.g., species plasticity and the existence of ecotypes; Randin et al., 2006). We exclude environmental extrapolation as a reason for the on average low performance of the predictions of the global model for Egypt, as we included environmental data for the area of Egypt in these models (but not the records), and hence, the predictions are not outside the realm of the global model (and hence do not represent an extrapolation).
While it is advisable to check for collinearity at training and prediction scales (Elith, Kearney, & Phillips, 2010), it is not always easy to maintain a representative set of variables that are uncorrelated at both scales. Although we minimized the correlation between environmental variables at global and regional scales to avoid unnecessarily high variance in model parameters, the correlation among environmental variables is, inevitably, not constant over space (Dormann et al., 2013). Some of the variables used at the global scale have high correlation in Egypt, making the reliability of predictions in Egypt less stable (Dormann et al., 2013; Elith et al., 2010). Furthermore, the quality of environmental variables is not constant in space. For example, the WorldClim data (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005; the source of most of the environmental variables used in this study) were adroitly prepared using interpolation of data from global weather stations. Weather stations are not evenly distributed in space: Climate data for areas such as Arabia and the Sahara (including Egypt) are interpolated using very few weather stations with high spatial clustering (see figure 1 in Hijmans et al., 2005), and hence, the interpolations are of potentially higher uncertainty that can affect the quality of calibrated models (Phillips, Anderson, & Schapire, 2006). This problem is not exclusive to the WorldClim data, but holds for any environmental layers derived from spatially-biased weather stations.
The environmental variables used may have been insufficient to characterize the species niche (Phillips et al., 2006). It is recommended to use proximal predictors (e.g., food sources or suitable roosting sites for bats) that directly describe the required resources and physiological limits than more indirect distal predictors (e.g., altitude; Austin, 2007; Merow et al., 2014). The use of proximal variables increases the transferability of models in space (Elith & Leathwick, 2009; Franklin, 2009). However, determining a set of species-specific proximal predictors is not possible without detailed knowledge of the ecology and physiology of each species, either unknown for most species (especially for bats) or not yet available at large scales (e.g., abundance of prey; Merow et al., 2014; Herkt, Matthias, Barnikel, Skidmore, & Fahr, 2016; Petitpierre, Broennimann, Kueffer, Daehler, & Guisan, 2017). The majority of SDM studies use (the easier to obtain) distal variables as surrogates for proximal variables; however, even if distal variables can indirectly describe the species requirements, the correlation between proximal and distal variables is not constant in space (Dormann et al., 2013; Elith & Leathwick, 2009; Merow et al., 2014). Examples of missed variables which can potentially improve model transferability for bats include locations of suitable roosting and foraging sites, proximity to water, food sources (Herkt et al., 2016; Razgour, Rebelo, Di Febbraro, & Russo, 2016). Regional models were calibrated for a limited environmental range (Figure S1), potentially contributing to the disagreement between regional and global model predictions.
While excessive model complexity can lead to overfitting to training data and consequent limited model transferability in space and time, we reject overfitting as a reason for the limited usefulness of priors. We limited overfitting using regularized modeling approaches, calibrated by spatial cross-validation blocks in a way that balances the number of presence locations and environmental variability between cross-validation folds (avoiding extrapolation) and adequately constrains the complexity of (both regional and global) models. That said, it is not clear how much model complexity optimization is affected by the limited number and quality of records (including sampling bias).
Predictions from global models interpolated to Egypt may well still describe the potential distribution of bats in Egypt. Their limited usefulness in our study only shows that the global dimension did not add new information, given the limitations of the available data from Egypt. If unbiased occurrence data were available, global models may indeed predict well in Egypt. Moreover, available bat records in Egypt are few and spatially-biased toward easily accessible areas, with the majority collected from relatively old literature and museum specimens. Most are opportunistic data gathered with an unknown sampling strategy (see Appendix S1). Due to their nocturnal and elusive behaviour, high maneuverability, and the need for specialized bat detectors for effective recording, it is challenging to obtain high-quality records for bats in developing countries (Razgour et al., 2016). Information on their geographical distribution is very limited, making bats highly under-represented in SDM studies (Herkt et al., 2016; Razgour et al., 2016), and Egypt is no exception. Finally, sampling bias can strongly affect model quality (Phillips et al., 2009), and while we attempted to correct for sampling bias in our models, we cannot quantify the efficiency of bias correction without bias-free data for comparison (Phillips et al., 2009; Warton et al., 2013), unavailable in most presence-only studies, especially in developing countries. The results of this study call for improved, systematic sampling of species occurrences in regions where currently only biased and scarce data are available.
We have shown that the use of global bat data did not improve regional model performance for Egypt. We relate this to the difficulty of calibrating SDMs of widespread species at extremely large study areas that cover many biogeographical regions and to data quality issues (mainly the quantity of available data dominated by high sampling bias). Due to the lack of high-quality data and limited environmental gradients in Egypt, regional SDMs seem to be insufficient to determine new survey sites (a point also made by Sánchez-Fernández et al., 2011). Improving the sampling of fauna and flora species from data-poor countries (such as Egypt, particularly from the less visited areas) would enhance regional SDMs in these countries and consolidate the usefulness of these models to discover new populations.
Although our results showed that predictions from global SDMs failed to improve regional predictions calibrated with low-quality and spatially-biased data, we still believe in great potential for SDMs that integrates global and regional data to improve future local sampling in data-poor countries like Egypt. Patterns of potential distribution (of global models interpolated to Egypt) can guide future surveys and help to discover new populations. In our analyses, we excluded Egyptian data for creating the global models to maintain consistency of comparisons between predictions of regional and global models. However, this is not necessary for real applications, and it would seem preferable to include regional data in a comprehensive model that covers the biogeographical region to improve model predictability. For example, to improve sampling of under-reported bat species in Egypt, we think that a larger-scale model should be created, with the study area determined objectively based on the available data from Egypt and adjacent arid areas (e.g., Arabia and the Sahara) in order to meet the stationarity assumption (constant species–environment relationships with no change in niche characteristics; Anderson & Gonzalez, 2011; Dormann et al., 2012) and then crop the prediction maps to Egypt. This is of mutual benefit not only for Egypt, but also for targetting efforts in the adjacent areas as well, which can help to improve the conservation status of some species. However, obtaining enough data from adjacent areas will remain challenging for many species.
ACKNOWLEDGMENTS We would like to express our sincere thanks to Petr Benda for comments on the global distribution of the bat species. An earlier version of the manuscript was improved by comments of Francis Gilbert and David R. Roberts. AE-G is sponsored by the German Academic Exchange Service (DAAD) through a GERLS scholarship. This work was partially performed on the computational resource “bwUniCluster” funded by the Ministry of Science, Research and Arts and the Universities of the State of Baden-Württemberg, Germany, within the framework program bwHPC. The article processing charge was funded by the German Research Foundation (DFG) and the University of Freiburg in the funding programme Open Access Publishing.
AUTHOR CONTRIBUTIONS AE-G and CFD contributed to idea and design of study, and comments and revisions; AE-G contributed to data curation and statistical analysis, and first drafted the writing. Both authors contributed critically to the drafts and gave final approval for publication.
CONFLICT OF INTEREST None declared.
Figure 1 The distribution of Asellia tridens at spatial (a) and environmental (b) space. The map a shows the species-specific global extent of this species, with dots representing the spatial distribution at global (blue) and regional (black) scales. Panel b shows a scatterplot of the first two PCA axes of all available environmental covariates within the entire study area. The first two axes account for 94.2% of the environmental variation. Blue and black dots are presence locations of the species outside and inside Egypt, respectively; light gray points are pixels without any sightings at global scale; dark gray points represent the available environmental space in Egypt. Figure S1 shows equivalent plot for all study species together
Figure 2 Mean permutation importance of environmental variables used at global (left) and regional (right) models (from Maxent). Dots and error bars represent the overall mean and standard deviation of the average permutation importance of the seventeen study species, respectively. Blue dots/bars represent environment-only models; red dots/bars represent comparable models with accessibility bias variables incorporated as predictors. When included, bias predictors have a high contribution (particularly distance to main cities at both scales, and distance to roads in Egypt), compared to many environmental variables. For more details on the environmental variables used, see Table S2
Figure 3 Boxplots for the geographical congruence (Schoener's D) between mean predictions of global and regional models for Egypt (with no priors). Schoener's D ranges from zero to one, representing situations of no to full congruence, respectively. “Env-only” are models calibrated only with environmental variables. “Bias-predictor” models add accessibility bias variables as predictors to the model. “Bias-corrected” models also use bias variables to set bias to zero during prediction (i.e., bias factored-out)
Figure 4 Boxplots for the mean AUC values (on cross-validation) calculated for different options of modeling algorithms, bias manipulations, and priors. (a) A comparison between mean AUC values of no-prior regional models and equivalent models that use different options of priors (without regional bias incorporated as predictors). (b) Same as a, with regional bias variables included as predictors
Figure 5 Geographical congruence between the predictions of regional SDMs calibrated without priors and the three versions of regional models that used a prior variable. Bias variables were not incorporated as predictors in the regional SDMs. There were three options of prior options: “Env-only” are predictions of global SDMs without incorporating sampling bias; “Bias-predictor” priors incorporate global accessibility bias variables as predictors in the model; and “Bias-corrected” priors incorporate bias-corrected (set to zero) predictions from global models for Egypt

Patrick Brading; Ahmed El-Gabbas; Samy Zalat; Francis Gilbert (2009). Biodiversity economics: the value of pollination services to Egypt. Egyptian Journal of Biology, 11, 45-51. Ahmed El-Gabbas, Carsten F. Dormann (2018). Wrong, but useful: regional species distribution models may not be improved by range-wide data under biased sampling. Ecology & Evolution 8(4):2196-2206. Ahmed El-Gabbas; Carsten F. Dormann (2018). Improved species-occurrence predictions in data-poor regions: using large-scale data and bias correction with down-weighted Poisson regression and Maxent. Ecography 41(7): 1161-1172. Ahmed El-Gabbas; Sherif Baha El Din; Samy Zalat; Francis Gilbert (2016). Conserving Egypt's reptiles under climate change. Journal of Arid Environments, 127:211-221. Ahmed El-Gabbas; Francis Gilbert (2016). The Desert Beauty Calopieris eulimene: a butterfly new to Egypt (Insecta: Lepidoptera). Zoology in the Middle, East 62:3(279-281). Tim Newbold; Tom Reader; Samy Zalat; Ahmed El-Gabbas, Francis Gilbert (2009). Effect of characteristics of butterfly species on the accuracy of distribution models in an arid environment. Biodiversity & Conservation 18: 3629-41. Tim Newbold; Francis Gilbert; Samy Zalat; Ahmed El-Gabbas, Tom Reader (2009). Climate-based models of spatial patterns of species richness in Egypt's butterfly and mammal fauna. Journal of Biogeography 36: 2085-95. Tim Newbold; Tom Reader; Ahmed El-Gabbas; Wiebke Berg; Wael M. Shohdi; Samy Zalat; Sherif Baha El Din; Francis Gilbert (2010). Testing the accuracy of species distribution models using species records from a new field survey. Oikos, 119(8): 1326-1334.