Modelling the spatial distribution of snake species in northwestern Tunisia using maximum entropy(Maxent)and Geographic Information System(GIS)

Mohsen Kalboussi·Hammadi Achour,2

Introduction

Snakes are among the least known reptile species in Tunisia.Little is known of their diversity and distribution(Domergue 1959;Nouira et al.1995),and some species have not been recorded for more than 50 years(Brito et al.2008).These animals are often killed when found in the field,even in some protected areas and people do not consider them as animals of special interest for conservation.Kroumiria has recorded snakes of six species(Nouira et al.1995)and three families,namely Natricidae,Colubridae and Lamprophiidae.Natricidae,is represented by two species:Natrix maura and Natrix astreptophora.Both of them inhabit watercourses and wetlands.N.maura has a wide geographical range in Tunisia(from the north to the oases in the south of the country;Nouira et al.1995).N.astreptophora,is a rare species,and is restricted to Kroumiria and its surrounding areas,after the compression of its range.The second family,Colubridae,is represented by three species:Hemorrhois hippocrepis,Coronella girondica and Macroprotodon mauritanicus.The horseshoe whip snake H.hippocrepis is known to occur in northern and central Tunisia.This species is a generalist in terms of habitat requirements(Segura et al.2007).C.girondica is a nocturnal species actually restricted to mountainous zones in northwestern Tunisia(Bons and Geniez 1996;Brito et al.2008).Its range seems to have contracted because no recent records have been made in other areas where it was known(Cap Bon peninsula;Domergue 1959).The false smooth snake M.mauritanicus is a habitat generalist and is widely distributed across Tunisia(Nouira et al.1995).It is also nocturnal.The Lamprophiidae are represented by Malpolon insignitus.It is the largest snake species in the region,inhabits a variety of environments,and has a wide distribution in the country(Nouira et al.1995).The northwestern part of Tunisia harbours some snake species at the limit of their distribution ranges in northern Africa(Bons and Geniez,1996;Schleich et al.1996),and some of them are rarely observed or reported(Brito et al.2008;Joger 2003;Nouira et al.1995).The ecology,habitat requirements and geographical distribution of many of these species are poorly known,mainly for those con fined to northwestern Tunisia.The scarcity of certain species needs urgent conservation efforts,beginning with greater understanding of their ecological requirements.

Species distribution modelling(SDM)aims at predicting potentialsuitablehabitatsacrossalandscapeusingempirical relationships between species distribution and environmental variables(Soberón and Peterson 2005;Thuiller and Münkemüller 2010).The recent development of techniques,including remote sensing,Geographic Information System(GIS),and data mining combined with increasing availability of large-scale environmental information in digital format have promoted signi ficant improvements in SDM techniques(Boyd and Foody 2011).Sillero and Tarroso(2010)listed several websites providing the most important free data for spatial analysis in ecology and provide also a list of the most commonly used functions in GIS analysis.SDM can be performed using various techniques and algorithms(Carpenter et al.1993;Busby 1991;Hirzel et al.2002;Stockwell and Peters 1999;Garcia-Rosello et al.2013;de Souza Muñoz et al.2011;Lehmann et al.2003;Thuiller 2003;Thuiller et al.2009;Leathwick et al.2006;Pearson et al.2002).Table 1 lists the most used predictive SDM packages.They can be classi fied and compared considering various aspects(i.e.,Presence-absence vs.presence-only modelling methods;Statistical vs.Machine learning modelling methods;speciesvs.community approach;Guisan and Zimmermann 2000;Zaniewski et al.2002).Each technique has its exclusive advantages and disadvantages.For a comprehensive review,see Brotons et al.(2004)and Elith et al.(2006).

In the present study,the maximum entropy(Maxent;Phillips et al.2006)model was used to predict the potential distribution range of our focal species for various reasons.First,the model has been successfully used for previous animal habitat modelling applications with presence-only data(Anadón et al.2012;Suarez-Seoane et al.2008;Wen et al.2015).Second,it has been speci fically recommended for sparse data con firming presence(Gibson et al.2007).According to Kaliontzopoulou et al.(2008)and Wisz et al.(2008),Maxent is quite robust to variations in sample size and produces useful results with sample sizes as small as 30 occurrences.Moreover,accordingtothe findingsofGraham et al.(2008),SDM techniques are not equally in fluenced by positional error in organism occurrence data.These authors carriedouttheirstudywithacellsizeof1 km2and0.01 km2for some regions,and incorporated locality spatial uncertainty(positional error)by using a random 5 km data shift radius.They concluded that certain methods,such as Maxent andboostedregressiontree(BRT)techniques,are robust to positional error.Finally,the model was designed to integrate with GIS software,thus making data easy to handle.However,sampling bias should be taken into account when using SDM techniques and particularly Maxent for modelling purposes.According to Yackulic et al.(2013),87%of Maxent models were based on data suffering from sample selection bias.Sampling bias can be addressed by reducing the number of occurrence records in oversampled regions using spatial filtering(Veloz 2009)or by using background data instead of using presence points data(Syfert et al.2013).Kramer-Schadt et al.2013 investigated the ef ficiency of spatial filtering versus background manipulation to reduce overprediction or underprediction in speci fic areas.They found that spatial filtering minimized both omission and commission errors.When sample size is insuf ficient to allow spatial filtering,they recommended the use of background data.

The present study is the first undertaken at a local scale in Kroumiria of our focal species.If other studies exist,they address much larger areas than this region,even at meso or global scales(Brito et al.2011).Local scale is used here to characterize areas equal to or smaller than 10,000 km2(Guisan and Hofer 2003),which corresponds to geographical scales in the order of 1/100,000 or less.However,the restricted range of species such as N.astreptophora and C.girondica,which are limited to northwestern Tunisia at the most southeastern limit of their global range(Schleich et al.1996),makes them particularly vulnerable to habitat degradationandclimatechange.Ouraimswere(1)toupdate the geographical distribution of the studied snake species in Kroumiriaand(2)toidentifythepredictorvariablesthatbest explained their distribution in northwestern Tunisia,by using presence data only.

Materials and methods

Study area

Fig.1 Digital elevation model and location of the study area showing presence data of the different snake families.The administrative limits within Kroumiria figure inside the map

Located in the northwest of Tunisia,the Kroumiria region is mountainous with extensive forests of cork oak(Quercus suber), and is the wettest region in Tunisia(1000–1500 mm peryear).Itcovers approximately 1775 km2.It extends south of the Mediterranean Sea and north of Wadi Medjerda and east from the Algerian border to Jebel El-Abiod(Fig.1).Topography is characterized by a marked elevation gradient ranging from 0-1200 masl with heterogeneous geology dominated by clayey sandstone rocks.Terrain gradient,computed from a Digital Elevation Model(DEM),ranges from 0°to 49°with an average of 17° and standard deviation of 10°.The climate is Mediterranean with annual mean temperature of 18.2°C(1975–2004).Maximum and minimum temperatures are 34.4 °C and 5.6 °C,respectively.Average annual rainfall is 912 mm,with 77%of total precipitation occurring in autumn and winter and only 4%in summer(Mechergui et al.2013).These climatic and physiographic conditions have favored the development of important and rich natural vegetation.Nowadays,the forest cover is dominated at low elevations(＜100 m)by kermes oak(Quercus coccifera),exceeding a height of eight meters in some areas.Once degraded,the kermes oak forest is replaced by a mosaic of natural vegetation dominated by Juniperus phoenicea,Retama raetam and Oleo-Lentiscetum associations.At medium elevations(300–800 masl),the forests mostly consist of pure stands of cork oak,and at higher elevations(＞800 m),by zeen oak(Quercus faginea)associations.In the north-western border of the region,the natural vegetation is dominated by pure maritime pine(Pinus pinaster)forests(Posner 1988).Scattered within the region,many sectors were planted by Pinus pinea,Eucalyptus sp.and maritime pine.Besides the preceding associations,we can also add some speci fic ecosystems,mainly perennial meadows and wetlands.Many dams were built in the region,in order to provide drinking water for Tunisian cities(located out of the Kroumiria)and for irrigation.

Species occurrence records

Species point locality data were obtained from intensive surveys conducted in the Kroumiria during 2000–2015.The field effort corresponded to 129 days,corresponding at least to 762 h.The number of individuals sighted by species are listed in Table 2.The most abundant species in the region were N.maura and M.mauritanicus,while C.girondica and N.astreptophora were least common.To our knowledge,these data present the most extensive records of the latter two species in Tunisia.

The reports included records of the latitude and longitude coordinates determined using Global PositioningSystem(GPS)with a positional accuracy between 10 m and 20 m.Records with multiple observations of the same species at the same location and time were treated as single observations.The georeferenced species point occurrence contained 63 observations.The difference between the number of individuals seen in the field and the number of species point locations is related to the fact that at some sites,more than one individual was found.We surveyed throughout the year,during times of snake activity or hibernation.During hibernation we searched under stones,piles of dead wood,trees bark,and fallen logs.During field surveys we identi fied,photographed,and released all living snakes,and identi fied and photographed dead snakes where possible.We recorded geographic coordinates at snake locations and brie fly described habitats and microhabitats.We searched for nocturnal snakes during the day in their retreats(no nocturnal field trips were conducted).As reptiles are ectotherms,surveys covered open areas(trenches,clearings,forests without undergrowth).Small archaeological sites scattered within forests,and some speci fic zones were avoided(mainly farmland and surroundings of human settlements).Almost all the natural ecosystems were visited more than once.Non-visited areas consisted mainly of inaccessible forests or border zones.

Table 2 Number of individuals seen in the field,by species

Environmental layers preparation

We used a set of environmental variables based on the available data,knowledge of species ecology and factors affecting species distribution(Anadón et al.2012;Bombi et al.2009;Brito et al.2011;Sow et al.2014).We listed all datasetsusedinthisstudyandtheircharacteristics(Table 3).To avoid including highly correlated variables in the model,we screened all variables for pairwise correlation using Pearson’s correlation coef ficient(r)for all candidate variables.To do this,we used the tool‘Raster Correlations and Summary Statistics’implemented in SDMToolbox(Brown 2014).Weconsideredthe variables tobehighly correlatedif r＞0.75 or r＜-0.75(Wen et al.2015).The results were considered signi ficant at p＜0.05.The Pearson correlation matrix did not show any signi ficant correlations,with only a few covariates exceeding a correlation coef ficient of 0.5(Table 4).In summary,the method of variable selectionresulted in seven variables for potential use in modelling.These include the following thematic layers:mean annual precipitation,elevation,slope gradient,aspect,distance to watercourses,Land Surface Temperature(LST)and Normalized Differential Vegetation Index(NDVI).

Table 3 List and characteristics of data used in this study

Daily precipitation data were obtained from government’s General Directorate of Water Resources for the period 1980–2010 for 49 weather stations.To increase the climatic data density and robustness of the analyses,some stations outside the Kroumiria region were used to interpolate precipitation data across the boundaries of the study area.After calculating the mean annual values,the precipitation data were imported to ArcGIS software for the exploratory spatial analysis.The ordinary kriging method was used to obtain estimates of average annual precipitation over a 30 year period.

Elevation data were generated from the 1 arc-second DEM(approximately 30 m)derived from the recently and freely available Shuttle RadarTopography Mission(SRTM)and retrieved from the United States Geological Survey(USGS)website(http://www.earthexplorer.usgs.gov). To cover the entire study area, the tile‘n36_e008_1arc_v3’was downloaded and then pit- filled using ArcGIS generic tools(ESRI,Inc.,Redlands,CA,USA).From this DEM,we computed slope gradient and aspect grids using a spatial analysis tool.Furthermore,as hydrology plays an important role in shaping the habitat of species,we also computed,using Euclidean distance tool,the thematic layer ‘distance to watercourses’from the stream lines digitized from topographic maps at the scale of 1:50,000.

Table 4 Correlation matrix of the selected variables

Images derived from the Landsat 8 Operational Land Imager(OLI)/Thermal Infrared Sensor(TIRS)satellite system(freely available at http://www.earthexplorer.usgs.gov)were used to derive Land Surface Temperature(LST)and Normalized Differential Vegetation Index(NDVI).The OLI sensor collects nine shortwave spectral bands over a 190 km swath with a 30 m spatial resolution,except the 15 m panchromatic band.The thermal bands 10(10.60–11.19 micrometers)and 11(11.50–12.51 micrometers),captured by the TIRS sensor,are initially collected at 100 meters and then resampled to 30 meters to match OLI multispectral bands(USGS 2015).To compute LST values,we followed all the steps of the process outlined in Landsat 8(L8)Data Users Handbook(USGS 2015).Brie fly,the approach involves four main calculation steps(Senay et al.2016):(1)Digital number values conversion to spectral radiance(2)spectral radiance conversion to brightness temperature;(3)surface emissivity computation;and(4)land surface temperature computation.All the above steps were computed using Map algebra tool.

Maximum entropy implementation

Maxent software version 3.3.3e,based on a maximum entropy algorithm(Phillips et al.2006),was used to prepare a statistical predictive model for the potential distribution of studied species.This model takes as input a set of layers or environmental variables,as well as a set of occurrence locations,and produces a model of the range of the given species.Maxent requires that all environmental data be in ASCII grids of the same resolution,extent and pixel alignment.As stated earlier,the input environmental data had spatial resolutions ranging from 30 to 60 m.Having the model developed at the resolution of the coarsest input was not practical for a species with local scale habitat requirements.Therefore a spatial resolution of 30 meters was chosen.The SDMToolbox for ArcGIS was used to process all raster inputs for use in Maxent.

As suggested by Pearson(2010),application of a model will have a little merit if the accuracy of the prediction is not assessed.The prediction results from Maxent were evaluated by a threshold-independent receiver operating characteristic(ROC)area under curve(AUC)values,calculated within the program.ROC curves plot true-positive rate against false-positive rate(Phillips et al.2006)and the area under the curve(AUC)was used as a measure of the overall fit of the model.AUC values＜0.5 are equal to those in a random prediction model;values＜0.7 indicate poor prediction ability,values 0.7–0.9 indicate moderate prediction ability;and values＞0.9 indicate strong predictive ability(Hemsing 2010).To improve prediction,the Maxent analysis was repeated 1000 times.Among the three replicate run type options(Subsample,Crossvalidate,and Bootstrap)available in Maxent,the ‘Subsample’method was used to sub-sample the presence data sets to evaluate model performance.Distribution maps were generated from the averages of the 1000 runs.

Results and discussion

Figure 2 shows the ROC and AUC values for our six studied species.For H.hippocrepis and M.insignitus,AUC values were 0.884 and 0.870,respectively,indicating moderate prediction ability.However,the AUC values for the other four studied species were＞0.9.Note however,that there was no correlation between AUC values and the number of presence points for all the studied species(r=0.318).Indeed,three species had the same number of presence points(C.girondica,H.hippocrepis and M.insignitus),but different AUC values(0.940,0.884 and 0.870,respectively).Explanation of this finding may be related to the spatial distribution of the studied species.High AUC values tend to result when species distribution is continuous.In contrast,AUC values tend to be low when species distribution is discontinuous.The frequency of presence data can also in fluence AUC values.These facts corroborate the prediction of Lobo et al.(2008),who concluded that the spatial distributions of rare species are usually better predicted than are the distributions of widespread species(case of C.girondica and N.astreptophora).

Fig.2 ROC curves for the considered snake species:a Coronella girondica,b Hemorrhois hippocrepis,c Macroprotodon mauritanicus,d Malpolon insignitus,e Natrix maura,f Natrix astreptophora

Fig.3 Jackknife tests for the six snake species:a Coronella girondica,b Hemorrhois hippocrepis,c Macroprotodon mauritanicus,d Malpolon insignitus,e Natrix maura,f Natrix astreptophora

Table 5 Relative contributions(%)and permutation importance(%)of the environmental variables to the Maxent model

Figure 3 and Table 5 show the results of jackknife tests.The most important factor explaining the presence of Ophidians was ‘distance to streams’,with ＞60%contribution,followed by elevation(＞10%contribution).Land surface temperature strongly has a relatively signi ficant effect in fluenced the spatial distribution only of C.girondica.The remaining variables had minor effects.Aspect and slope made the least contribution.Our findings did not corroborate those of Anadón et al.(2007,2012)or Santos et al.(2009),who found that climatic variables were the most important factors in fluencing the presence of their studied reptile species.In his case study of four snake species in Spain,Muthoni(2010)found that climatic variables,especially temperature-related gradients,predicted signi ficantly better than biophysical variables.Silva-Rocha et al.(2015),in their study on invasive species in the Balearic Islands,concluded thatclimatic variables,including precipitation of the warmest quarter and seasonality of precipitation,had a determinant in fluence on the distribution of M.mauritanicus.A possible explanation of these findings is that biophysical variables were not included in our modelling.According to Guisan and Hofer(2003),biophysical variables at local scale better explain the spatial distribution of species than climatic variables,which are explanatory at meso-scale.This may be attributed to the fact that climate is relatively stable over local areas such as Kroumiria,thus our generated model was more strongly in fluenced by biophysical variables.Kroumiria is the wettest region of Tunisia and is drained by many watercourses with seasonal flow;it is then obvious that snakes prefer such habitats because of their high biological productivity.The distributions of N.astreptophora and N.maura are linked to the presence of water(Harris et al.2008;Rouag and Benyacoub 2006).However,H.hippocrepis and M.insignitus are generalist species in terms of habitat requirements(Segura et al.2007).They were found near watercourses,simply because these habitats provided suitable conditions for both of them(presence of shelters,food resources).

Figure 4 shows the predicted distribution areas of the six studied snake species.The habitat-suitability map of C.girondica shows a probability of presence ranging from 0(unsuitable habitat)to 0.91(optimal habitat)with an average of 0.33 and a standard deviation of 0.19.Areas of high habitat suitability dominated mostly the northern and southwestern parts of the study area(El Feidja National Park and El Ghorra mountain),forming two large patches,which correspond mainly to the mountainous and elevated locations.Wide regions in both eastern and southeastern parts of Kroumiria appeared almost unsuitable as habitats for the species.The total extent of suitable habitat(probability≥0.7)was approximately 83 km2,which represents 4.6%of the entire study area.

Coronella girondica,while restricted to northwestern Tunisia is not limited to Kroumiria.Its range is limited to holm and zeen oak forests which extend out of the study area to its continuous southwestern regions(Nefza)where we found it.It is rare in this part of the country.The elevation range of C.girondica was between 47 m and 1020 masl.The last finding contradicts what is known about the ecology of the species(an inhabitant of mountainous areas;Brito et al.2008).Only Pastorelli et al.(2007)reported a specimen caught in Italy at low elevation(47 masl)in 1855,while Schleich et al.(1996)stated that C.girondica occupies Mediterranean lowlands but did not specify localities or elevations.The species seems to live in a wide range of elevations,but tends to occupy elevated areas because of the transformations of its natural environment,since it is known to live in unaltered habitats(Segura et al.2007).

The potential distribution model of H.hippocrepis showed a probability of presence ranging from 0.12 to 0.76,with an average of 0.49.The potentially suitable areas were mainly concentrated in the northern part of the study area.The geographical extent of suitable habitat(probability≥0.7)was about 71 km2,or 3.9%of Kroumiria.These patches correspond to the areas where the species has been recorded.Wide regions in both northern and southern Kroumiria appear completely unsuitable as habitats for H.hippocrepis despite its generalist character(Segura et al.2007).Our observations indicate that the species lives not only in natural environments(oak and pine forests,meadows,trenches),but also close to human settlements(abandoned buildings)and in olive plantations,and archaeological sites.It occupies elevations ranging from 27 to 629 masl.

The habitat suitability map of M.insignitus showed a probability of presence ranging from 0 to 0.77 with an average of 0.49 and a standard deviation of 0.12.The apparently suitable habitats were highly fragmented and discontinuous with areas of high suitability mostly in the northwestern part of the study area.The geographical extent of suitable habitat(probability≥0.7)was about 23 km2,or 2%of the study area.

The apparently suitable habitats of M.mauritanicus appeared to be highly fragmented,with areas of high suitability mostly located in the western part of the region,but also in the northern mountain chains.The total extent of apparently suitable habitat(probability≥0.7)was 72.8 km2,or 4%of the study area.Overall,seven patches were identi fied throughout the region.Among these,one large patch(＞26 km2)was clearly visible in the generated map and corresponded to El Feidja National Park and its surrounding areas.The species reaches the highest elevation of all Ophidians(1070 masl.).

The habitat suitability of N.astreptophora showed a probability of presence ranging from 0 to 0.89,with an average of 0.44 and a standard deviation of 0.16.Its geographical extent was about 68 km2(3.8%of the study area),in the northern part of the Kroumiria.Areas of high habitat suitability,although without records of the species,were observed mostly in both northeastern and southeastern parts of Kroumiria,suggesting the power of the model to detect suitable areas in unsurveyed regions.Our work is the first to deal with this species after its taxonomic separation from its congener,Natrix natrix(Pokrant et al.2016).

Fig.4 Predicted distribution areas for the studied species

In Tunisia,N.astreptophora is restricted to the northwestern part of the country(Kroumiria and adjacent areas,i.e.Nefza).We found it once near human settlements.Its main habitats correspond to forests bordering meadows,cork and zeen oak forests and Eucalyptus plantations.It is active not only during warm seasons,but also during winter on hot days.The elevation range of N.astreptophora ranges from a few meters above sea level to 1002 masl on El Ghorra mountain.

The potential distribution of N.maura appeared to be more restricted than the other species and it was predicted to be concentrated in the northern part of Kroumiria.The total extent of apparently suitable habitats(probability≥0.7)was 40.5 km2,or 2.3%of the study area.The model of suitable habitat suggested that despite the obviously wide ecological tolerance of N.maura,there were areas of unsuitable habitats within its natural distribution.According to our observations,N.maura can have common hibernacula and it was common to find several individuals sharing the same shelter or under wintering retreats very close to one another.Notice,however,that no vipers were encountered during any field survey.We searched in particular for Vipera latastei in the areas where it was reported to occur(Chpakowsky and Chnéour 1953).Based on this result,we follow Brito et al.(2008)in concluding that the species is now extinct in Tunisia.

Conclusions

Our results showed that ‘distance from streams’was the main factor that set the distributional limits of these species at fine spatial scale.However,it is probable that other climatic and non-climatic factors(biotic interactions,dispersal abilities,and evolutionary adaptations)that were not incorporated in our models in fluence the spatial distribution of these species.Among the studied snakes,N.astreptophora and C.girondica are the main species of conservation interest and could be considered as focal species requiring priority conservation measures.Both species are mainly located in the Kroumiria region and have restricted habitatrequirements.Assome suitable habitats for both species correspond to small areas in the region,their conservation may reinforce local populations and,then,enhance the preservation of both of them.

Acknowledgements We thank all the persons,mainly students,who helped during field trips.We also thank reviewers for their constructive comments to improve this paper.

Anadón JD,Giménez A,Martínez M,Palazón JA,Esteve MA(2007)Assessing changes in habitat quality due to land use changes in the spur-thighed tortoise Testudo graeca using hierarchical predictive habitat models.Divers Distrib 13:324–331

Anadón JD,Giménez A,GraciáE,Perez I,Ferrández M,Fahd S,El Mouden H,Kalboussi M,Jdeidi T,Larbes S,Rouag R,Slimani T,Znari M,Fritz U(2012)Distribution of Testudo graeca in the western Mediterranean according to climatic factors.Amphib Reptil 33:285–296

Bombi P,Luiselli L,Capula M,Salvi D(2009)Predicting elusiveness:potential distribution model of the southern smooth snake,Coronella girondica,in Italy.Acta Herpetol 4:7–13

Bons J,Geniez P(1996)Amphibiens et Reptiles du Maroc(Sahara Occidental compris).Atlas biogéographique.Associación Herpetológica Española,Barcelona

Boyd D,Foody GM(2011)An overview of recent remote sensing and GIS based research in ecological informatics.Ecol Inform 6:25–36

Brito JC,Feriche M,Herrera T,Kaliontzopoulou A,Martínez-Freiría F,Nesbitt D,Omolo D,Ontiveros D,Quiñoz L,Pleguezuelos JM,Santos X,Sillero N(2008)En los límites de su distribución:an fibios y reptiles paleárticos en el noroeste de Tu´nez.Bol Asoc Herpetol Esp 19:75–82

Brito JC,Fahd S,Geniez P,Martinez-Freiría F,Pleguezuelos JM,Trape JF(2011)Biogeography and conservation of viperids from North-West Africa:an application of ecological niche-based models and GIS.J Arid Environ 75:1029–1037

Brotons L,Thuiller W,Arau´jo MB,Hirzel AH(2004)Presenceabsence versus presence-only modelling methods for predicting bird habitat suitability.Ecography 27:437–448

Brown JL(2014)SDMtoolbox:a python-based GIS toolkit for landscape genetic,biogeographic and species distribution model analyses.Methods Ecol Evol 5:694–700

Busby JR(1991)BIOCLIM—a bioclimatic analysis and prediction system.In:Margules CR,Austin MP(eds)Nature conservation:cost effective biological surveys and data analysis.CISIRO,Canberra,pp 64–68

Carpenter G,Gillison AN,Winter J(1993)DOMAIN:a flexible modelling procedure for mapping potential distributions of plants and animals.Biodivers Conserv 2:667–680

Chpakowsky N,Chnéour A(1953)Les serpents de Tunisie.Bull Soc Sci Nat Tunis 6:125–146

de Souza Muñoz ME,De Giovanni R,de Siqueira MF,Sutton T,Brewer P,Pereira RS,Canhos DAL,Canhos VP(2011)openModeller:a generic approach to species’potential distribution modelling.GeoInformatica 15:111–135

Domergue CA(1959)Liste des ophidiens de Tunisie,de l’Algérie et du Maroc.Arch Inst Pasteur Tunis 36:157–161

Elith J,Graham C,Anderson R,Dudik M,Ferrier S,Guisan A,Hijmans R,Huettmann F,Leathwick J,Lehmann A,Li J,Lohmann L,Loiselle B,Manion G,Moritz C,Nakamura M,Nakazawa Y,Overton J,Peterson A,Phillips S,Richardson K,Scachetti-Pereira R,Schapire R,Soberon J,Williams S,Wisz M,Zimmermann N(2006)Novel methods improve prediction of species’distributions from occurrence data. Ecography 29:129–151

Garcia-Rosello E,Guisande C,Gonzalez-Dacosta J,Heine J,Pelayo-Villamil P,Manjarras-Hernandez A,Vaamonde A,Granado-Lorencio C(2013)ModestR:a software tool for managing and analyzing speciesdistribution map databases.Ecography 36:1202–1207

Gibson L,Barrett B,Burbidge A(2007)Dealing with uncertain absences in habitat modelling:a case study of a rare grounddwelling parrot.Divers Distrib 13:704–713

Graham CH,Elith J,Hijmans RJ,Guisan A,Townsend Peterson A,Loiselle BA,Anderson RP,Dudk M,Ferrier S,Huettmann F,Leathwick J,Lehmann A,Li J,Lohmann L,Loiselle B,Manion G,Moritz C,Nakamura M,Nakazawa Y,Overton J,Phillips S,Richardson K,Pereira RS,Schapire R,Soberón J,Williams S,Wisz M,Zimmermann N(2008)The in fluence of spatial errors in species occurrence data used in distribution models.J Appl Ecol 45:239–247

Guisan A,Hofer U(2003)Predicting reptile distributions at the mesoscale:relation to climate and topography.J Biogeogr 30:1233–1243

Guisan A,Zimmermann NE(2000)Predictive habitat distribution models in ecology.Ecol Model 135:147–186

Harris DJ,Carretero MA,Brito JC,Kaliontzopoulou A,Pinho C,Perera A,Vasconcelos R,Barata M,Barbosa D,Carvalho S,Fonseca MM,Pérez-Lanuza G,Rato C(2008)Data on the distribution of the terrestrial herpetofauna of Morocco:records from 2001–2006.Herpetol Bull 103:19–28

Hemsing LØ(2010)GIS modelling of potential natural vegetation(PNV):a methodological case study from south-central Norway.Master thesis,Norwegian University of Life Sciences.https://brage.bibsys.no.Accessed 02.09.16

Hirzel AH,Hausser J,Chessel D,Perrin N(2002)Ecological-niche factor analysis:how to compute habitat-suitability maps without absence data?Ecology 83:2027–2036

Joger U(2003)Reptiles and amphibians of southern Tunisia.Kaupia 12:71–88

Kaliontzopoulou A,Brito JC,Carretero MA,Larbes S,Harris DJ(2008)Modelling the partially unknown distribution of wall lizards(Podarcis)in North Africa:ecological af finities,potential areas of occurrence,and methodological constraints.Can J Zool 86:992–1001

Kramer-Schadt S,Niedballa J,Pilgrim JD,Schroder B,Lindenborn J,Reinfelder V,Stillfried M,Heckmann I,Scharf AK,Augeri DM,Cheyne SM,Hearn AJ,Ross J,Macdonald DW,Mathai J,Eaton J,Marshall AJ,Semiadi G,Rustam R,Bernard H,Alfred R,Samejima H,Duckworth JW,Breitenmoser-Wuersten C,Belant JL,Hofer H,Wilting A(2013)The importance of correcting for sampling bias in MaxEnt species distribution models.Divers Distrib 19:1366–1379

Leathwick JR,Elith J,Hastie T(2006)Comparative performance of generalized additive models and multivariate adaptive regression splines for statistical modelling of species distributions.Ecol Model 199:188–196

Lehmann A,Overton JM,Leathwick JR(2003)GRASP:generalized regression analysis and spatialprediction.EcolModel 160:165–183

Lobo JM,Jiménez-Valverde A,Real R(2008)AUC:a misleading measure of the performance of predictive distribution models.Global Ecol Biogeogr 17:145–151

Mechergui T,Pardos M,Boussaidi N,Hasnaoui B,Jacobs DF(2013)Development of cork oak(Quercus suber L.)seedlings in response to tree shelters and mulching in northwestern Tunisia.J For Res 24:193–204

Muthoni FK(2010)Modelling the spatial distribution of snake species under changing climate scenario in Spain.Master thesis,Faculty of Geo-information Science and Earth Observation

Nouira S,Blanc CP,Ktari MH(1995)Biodiversitéde l’herpétofaune tunisienne.I.Les ophidiens.Bull Soc Sci Nat Tunis 24:76–94

Pastorelli C,Laghi P,Melloni L(2007)Distribuzione di Coronella girondica(Daudin,1803)in Romagna(Reptilia:Squamata:Colubridae).Acta Herpetologica 2:121–217

Pearson RG(2010)Species’distribution modelling for conservation educators and practitioners.Lessons in conservation.American Museum of Natural History,pp 54–89.http://ncep.amnh.org.Accessed 05.09.16

Pearson RG,Dawson TP,Berry PM,Harrison PA(2002)SPECIES:a spatial evaluation of climate impact on the envelope of species.Ecol Model 154:289–300

Phillips SB,Aneja VP,Kang D,Arya SP(2006)Modelling and analysis of the atmospheric nitrogen deposition in North Carolina.Int J Glob Environ 6:231–252

Pokrant F,Kindler C,Ivanov M,Cheylan M,Geniez P,Böhme W,Fritz U(2016)Integrative taxonomy provides evidence for the species status of the Ibero-Maghrebian grass snake Natrix astreptophora.Biol J Lin Soc 118:873–888

Posner SD(1988)Biological diversity and tropical forests in Tunisia.Washington

Rouag R,Benyacoub S(2006)Inventaire et écologie des reptiles du Parc national d’El Kala(Algérie).Bull Soc Herpétol Fr 117:25–40

Santos X,Brito JC,Caro J,Abril AJ,Lorenzo M,Sillero N,Pleguezuelos JM(2009)Habitat suitability,threats and conservation of isolated populations of the smooth snake(Coronella austriaca)in the southern Iberian Peninsula.BiolCons 142:344–352

Schleich HH,Kästle W,Kabisch K(1996)Amphibians and Reptiles of North Africa.Koeltz Sci.ed.,Koenigstein

Segura C,Feriche M,Pleguezuelos JM,Santos X(2007)Specialist and generalist species in habitat use:implications for conservation assessment in snakes.J Nat Hist 41:2765–2774

Senay GB,Friedrichs M,Singh RK,Velpuri NM(2016)Evaluating Landsat 8 evapotranspiration for water use mapping in the Colorado River Basin.Remote Sens Environ.doi:10.1016/j.rse.2015.12.043

Sillero N,Tarroso P(2010)Free GIS for herpetologists:free data sources on Internet and comparison analysis of proprietary and free/open source software.Acta Herpetol 5:63–85

Silva-Rocha I,Salvi D,Sillero N,Mateo JA,Carretero MA(2015)Snakes on the Balearic Islands:an invasion tale with implications for native biodiversity conservation. PLoS ONE 10:e0121026.doi:10.1371/journal.pone.0121026

Soberón J,Peterson TA(2005)Interpretation of Models of Fundamental Ecological Niches and Species’Distributional Areas.Biodivers Inform 2:1–10

Sow AS,Martinez-Freiria F,Dieng H,Fahd S,Brito JC(2014)Biogeographical analysis of the Atlantic Sahara reptiles:environmental correlates of species distribution and vulnerability toclimate change.J Arid Environ 109:65–73

Stockwell D,Peters D(1999)The GARP modelling system:problems and solutions to automated spatial prediction.Int J Geogr Inf Sci 13:143–158

Suarez-Seoane S,Garcia de la Morena EL,Morales Prieto MB,Osborne PE,de Juana E(2008)Maximum entropy niche-based modelling of seasonal changes in little bustard(Tetrax tetrax)distribution.Ecol Model 219:17–29

Syfert MM,Smith MJ,Coomes DA(2013)The effects of sampling bias and model complexity on the predictive performance of MaxEnt species distribution models.PLoS ONE 8:e55158

Thuiller W(2003)BIOMOD—optimizing predictions of species distributions and projecting potential future shifts under global change.Glob Change Biol 9:1353–1362

Thuiller W,Münkemüller T(2010)Habitat suitability modelling.In:Møller AP,Fiedler W,Berthold P(eds)Effects of climate change on birds.Oxford University Press,New York,pp 77–85

Thuiller W,Lafourcade B,Engler R,Arau´jo MB(2009)BIOMOD—a platform for ensemble forecasting of species distributions.Ecography 32:369–373

USGS(2015)Landsat 8(L8)data users handbook.Earth Resources Observation and Science(EROS)Center 8,97

Veloz SD(2009)Spatially autocorrelated sampling falsely in flates measures of accuracy for presence-only niche models.J Biogeogr 36:2290–2299

Wen L,Saintilan N,Yang XH,Hunter S,Mawer D(2015)MODIS NDVI based metrics improve habitat suitability modelling in fragmented patchy floodplains.Remote Sens Appl Soc Environ 1:85–97

Wisz MS,Hijmans RJ,Li J,Peterson AT,Graham CH,Guisan A,Elith J,Dudík M,Ferrier S,Huettmann F,Leathwick JR,Lehmann A,Lohmann L,Loiselle BA,Manion G,Moritz C,Nakamura M,Nakazawa Y,Overton JM,Phillips SJ,Richardson KS,Scachetti-Pereira R,Schapire RE,Soberón J,Williams SE,Zimmermann NE(2008)Effects of sample size on the performance of species distribution models.Divers Distrib 14:763–773

Yackulic CB,Chandler R,Zipkin EF,Royle JA,Nichols JD,Campbell Grant EH,Veran S(2013)Presence-only modelling using MAXENT:when can we trust the inferences?Methods Ecol Evol 4:236–243

Zaniewski E,Lehmann A,Overton JM,McC J(2002)Predicting species spatial distributions using presence-only data:a case study of native New Zealand ferns.Ecol Model 157:261–280