Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Assessing the Spatial Scale Effect of Anthropogenic Factors on Species Distribution

Abstract

Patch context is a way to describe the effect that the surroundings exert on a landscape patch. Despite anthropogenic context alteration may affect species distributions by reducing the accessibility to suitable patches, species distribution modelling have rarely accounted for its effects explicitly. We propose a general framework to statistically detect the occurrence and the extent of such a factor, by combining presence-only data, spatial distribution models and information-theoretic model selection procedures. After having established the spatial resolution of the analysis on the basis of the species characteristics, a measure of anthropogenic alteration that can be quantified at increasing distance from each patch has to be defined. Then the distribution of the species is modelled under competing hypotheses: H0, assumes that the distribution is uninfluenced by the anthropogenic variables; H1, assumes the effect of alteration at the species scale (resolution); and H2, H3 … Hn add the effect of context alteration at increasing radii. Models are compared using the Akaike Information Criterion to establish the best hypothesis, and consequently the occurrence (if any) and the spatial scale of the anthropogenic effect. As a study case we analysed the distribution data of two insular lizards (one endemic and one naturalised) using four alternative hypotheses: no alteration (H0), alteration at the species scale (H1), alteration at two context scales (H2 and H3). H2 and H3 performed better than H0 and H1, highlighting the importance of context alteration. H2 performed better than H3, setting the spatial scale of the context at 1 km. The two species respond differently to context alteration, the introduced lizard being more tolerant than the endemic one. The proposed approach supplies reliably and interpretable results, uses easily available data on species distribution, and allows the assessing of the spatial scale at which human disturbance produces the heaviest effects.

Introduction

Human presence and activities inevitably produce changes in the surrounding environment that may affect the occurrence of a species in a given territory: some species may go locally extinct because they cannot find suitable conditions, while some others may spread to new territories because of the new conditions. Consequently, anthropogenic habitat alteration is considered as one of the dominant forces shaping the spatial distribution of species [1,2], with outcomes that are sometimes beneficial for the species and sometimes not. The conservation biology has grown and evolved in response to human threats to the natural world in order to set the best scientific guidance to preserve sites and species of special interest. One of the central topics for conservation biologists is to be able to reliably evaluate how the distribution of a species is affected by anthropogenic habitat alteration and loss. The issue is crucial and may be the starting point to elaborate really effective conservation plans.

In the last two decades the development of techniques which allow modelling the spatial distribution of a species on the basis of the relationship between species and environment has given a powerful tool to face the problem: an increasing number of studies have actually made use of spatial distribution models (SDM) – also known as ecological niche models [3] – to answer questions related to conservation tasks [46]. As their use has increased, greater and greater attention has been paid to the effects of the spatial scale (with different meanings) on the modelling process [714]. Besides the well-known importance of spatial scale in ecology [15,16], this particular attention comes from the plain consideration that the distribution of a species is the result of the combination of different factors (summarized in the BAM diagram – Figure 1 – where “B” stands for biotic interactions, “A” for abiotic conditions, and “M” for movement, i.e. area accessibility [3,17]), which act in the geographic space and are related each other in a scale-dependent way [11,18]. The majority of the studies focusing on spatial scale were devoted at analysing spatial grain and/or extent (sensu Wiens [15]) and how SDM respond to their changes [7,8,10,12,14]. These two components in the SDM studies have been generally associated respectively to the resolution (pixel size) and the geographic range (extent of the study area) at which the models are developed [8]. A relatively less studied issue, which is related to the “extent” side of spatial scale [15,19], was the relationship between the characteristics of the surroundings of a pixel on the predicted suitability of that pixel [9,15]. This aspect of extent (known as “patch context” in landscape ecology [19]) may play an important role in shaping the distribution of a species, for example by conditioning the dispersal ability [2022], i.e. the “M” factor of the BAM diagram. A support for the relevant role of patch context comes from the few studies that explicitly incorporate it among the variables employed to build the models [9,13,23]: all of them found a significant effect of patch context, and in some cases [9] they were able to assess the spatial extent to be considered “context”. Nevertheless, the experimental design adopted by all these studies was based on a sample of intensively surveyed sites, thus limiting the geographic range of the study area and preventing the assessment of general pattern of dependence between the species distribution and the spatial scale of the context. Furthermore, context was defined without an explicit reference to the anthropogenic alteration, thus none of the above studies explicitly addressed the link between human-induced alteration and species occurrence, nor the scale at which human activities act on species distribution. In the light of this lacks, our aim is to propose a generalizable approach to assess the occurrence and the spatial scale (in term of context extent) of the anthropogenic effect on the species distribution by combining species presence data, SDM and information-theoretic model comparison (IMTC) [24].

thumbnail
Figure 1. BAM diagram.

Simplified version of the Biotic-Abiotic-Movement diagram from Sauge et al. [16]. “G” represents the geographic space. “B” is the part of G that presents the correct set of biotic conditions: in this simplified version B is assumed not to constrain species distribution [3,16]. “A” is the part of G that holds suitable abiotic conditions. “M” represents the sub-area of G that has been accessible and explored by the species. The intersection between M and A defines the species distribution (G0). GI is the area that is potentially suitable, but has not been accessible to the species.

https://doi.org/10.1371/journal.pone.0067573.g001

SDM have been developed to predict the distribution of a species on the basis of known occurrences and a set of environmental layers (e.g., annual precipitation, annual mean temperature, land use): the model produces suitability scores for all the sites within a given area by capturing species-environment relationship [5]. Since links between species and environment are the key-step of the procedure, these models are also known as ecological niche models [3,6]. Without entering the debate about nomenclature and related concepts [3,2527], we adopt the term “SDM” according to Peterson and Soberón [3], since the objective of our modelling is the estimate of the actual distribution, taking into account accessibility (though indirectly; see below) [3]. Depending on the biological questions, SDM may be used as a predictive or explanatory tool: while in the former case SDM are employed to make predictions to new sites (either in space or in time) [28,29], in the latter case SDM are built to investigate the causes of the observed distribution, i.e., models serve to evaluate the effects of the variables used to build the SDM on the performance of the model itself: it is implicitly assumed that any effect on the model corresponds to a relationship between the variable and the distribution [30,31]. Our use of SDM belongs to this latter case. ITMC is a general well known approach in ecology and other branches of sciences [32]. It allows comparing different hypotheses, each represented by a model, searching for the best one [33]. The competing models are ranked on the basis of an information criterion (such as the Akaike Information Criterion or the Bayesian Information Criterion) which measures the amount of the information not captured by the model, weighed for the complexity of the model itself [32,33]. In addition to the ranking of the models from the best to the worst, it is also possible to evaluate how much a model, and its underlying causal hypothesis, is better than another one [33]. Our idea is to combine these two techniques in order to answer the initial question - how is the observed distribution of a species affected by anthropogenic habitat alteration? - paying attention to the spatial scale of the processes involved.

Methodological Scheme

The approach can be summarized as follows: i) establishing the resolution (grain component of the spatial scale) of the analysis, taking into consideration the characteristic of the species under study; ii) defining a measure of habitat alteration that can be applied at different spatial scales (i.e. quantifiable at increasing distance from each pixel); iii) obtaining a set of geographic layers each representing the alteration at a given scale iv) formulating a set of competing hypotheses with or without the habitat alteration factors; v) translating each hypothesis into an SDM; vi) comparing the obtained SDM using the ITMC, in order to assess the SDM (and the hypothesis underlying it) best explaining the species occurrence. In formulating the competing hypotheses, the null hypothesis assumes that anthropogenic changes do not affect distribution (H0): an approximation of H0 is a model which includes only climatic and topographical variables (that is to say a topoclimatic model), since it predicts species occurrence when habitat variables reflect their climatic potential and all the potentially suitable areas are accessible to the species [3]. Even if these assumptions may seem unrealistic, they serve as a starting point for comparing the models generated under hypotheses that call for anthropogenic effects. We are aware that other non-anthropogenic factors may actually affect both dispersal and pixel characteristics and consequently trim the potential distribution, but we consider that: i) a careful and expert-based choice of the extent of the study area may help producing a more realistic H0 [6,34]; ii) the aim is not to build the fittest SDM, but to assess if adding anthropogenic factors improves the fit of the model. We expect that anthropogenic alteration acts either on the quality of a pixel as on the permeability of its surroundings, thus reducing or modifying the potential distribution of a certain amount [35]. So we are interested in the trimming effect produced by an anthropogenic variable on the topoclimatic model.

Within the previous framework, each of the alternative HYPOTHESIS H1, H2 etc. should include the topoclimatic variables and one or more layers that describe anthropogenic alterations at different scales. In this way a set of nested models are obtained, each differing from H0 only for the presence of the human effect among variables. In order to illustrate our theoretical scheme at work and explain the detail of the procedure, we applied it to a set of data on two Sardinian lizard species.

Species and study area

We applied our scheme to the case of two congeneric lacertid lizards which inhabit a large Mediterranean island. The choice of a large island is motivated by two reasons: i) to bypass the problem of defining the limits of the study area; ii) to avoid the problem of resources limitation due to small dimension of the area.

Sardinia (40.10° N, 8.95° E, Italy) is a large island (24090 km2) located in the Tyrrhenian Sea (Western Mediterranean Sea). The numerous islets that surround the main island were not included, because of the particular conditions they represent. Following the Köppen-Geiger climate classification [36], the climate belongs to the Cs type (Warm temperate climate with dry summer), subtype a (hot summer) and b (warm summer). Vegetation consists of forests (primarily evergreen oak and cork trees and deciduous woods with oak and chestnut) and Mediterranean scrubs [37].

The island is inhabited by both endemic and introduced species [38], including the lacertid lizards of the genus Podarcis, the Tyrrhenian wall lizard P. tiliguerta (Gmelin, 1798), and the Italian wall lizard P. siculus (Rafinesque-Schmaltz, 1810). The former is a small lizard (snout-ventral length in general up to 6.5 cm [39,40]) endemic to Corsica and Sardinia and their surrounding islets. It shows a great ecological plasticity occurring preferentially in sparse scrub with rocky outcrops, but avoiding patches with too dense woodland and habitats characterized by intensive agriculture [40,41]. It ranges from the sea level up to about 1800 m a.s.l. [41,42]. The Italian wall lizard, is a medium-sized lizard, slightly larger than P. tiliguerta (snout-ventral length up to 9 cm [39]), occurring in southern-central Europe (continental Italy and on many islands, coastal region of Slovenia and Croatia and some areas of Montenegro). It shows a great colonization ability (naturalized populations are found in Portugal, Spain, France, Turkey, Tunisia, Libya and USA [39,43]), and it has been introduced in Sardinia in subsequent historical or protohistorical times, maybe following commercial routes from Sicily [44,45]. The Italian wall lizard can use many different habitats (including anthropic ones) characterized by sunny and open areas [39,4143]. These two species have been chosen because of their different response to anthropization, the Italian wall lizard being more tolerant. So we should expect differences in the effects of anthropogenic habitat alterations between them.

A total of 934 species occurrence points were collected during surveys from 2000 and 2009 as part of the study carried out by the Museo di Storia Naturale dell’ Università di Firenze, Sezione di Zoologia and the Sezione Sardegna of the Italian Society of Herpetology S.H.I., (Societas Herpetologica Italica) in order to produce the Italian Herpetological Atlas and Fauna d’Italia Reptilia volume [38,43] and by the Regione Autonoma della Sardegna for the realisation of the regional atlas (still on-going). This research was carried out under permits released by the Italian Ministry of Environment. None of the animals was captured and/or manipulated, and species identification was made by sight (this kind of study does not need permits under the current Italian and European legislation; furthermore, the sampled area falls outside any kind of restricted areas which needs permission). Geographical locations of species occurrence points were recorded using a GPS. No data about absence were available, so the sample represented a presence-only dataset. Data were resampled to a 100 × 100 m grid obtaining respectively 291 presence cells for P. tiliguerta and 355 for P. siculus. We chose this spatial grain on the basis of available data about the dispersal ability of the lizards of the genus Podarcis. In particular, maximum home range dimension has been demonstrated to be of about 300 m2 [4648] and homing ability has been registered up to 150 m distance [46,49]. So, a 100 m wide cell is expected to enclose the whole home range of a lizard and to represent its dispersal ability; for these reasons this scale can reliably represent the species spatial scale.

Two kinds of environmental data at the same spatial resolution (100 m) were directly obtained from the web: a digital elevation model (DEM) (http://srtm.csi.cgiar.org) and a map of the land use relating to 2008 (http://www.sardegnageoportale.it). Starting from this cartographic base, two sets of variables were generated: i) topographic, from DEM (altitude, slope and potential solar radiation), and ii) anthropic, from land use (see below). We had initially chosen not to use the many bioclimatic databases available on the web because none of them is detailed enough to allow 100 × 100 m resolution (e.g., the 19 bioclimatic variables of the Worldclim series, www.worldclim.org, has a maximum resolution of about 1 × 1 km at the latitude of Sardinia). Nevertheless, we compared the variables derived from DEM (resampled to 1 km) to those from Worldclim [50]: the high values of the correlation coefficients among most of them (12 out of 19; Table S1) supported the choice of using altitude, slope and potential solar radiation as surrogates of these bioclimatic variables. On the other hand, since some ecologically meaningful climatic variables were not well represented in our H0, the null model would have been easily improvable, via correlation, by any kind of appended variable. To avoid this issue, we decided to include also those bioclimatic variables that were not highly correlated with the topographic ones or among each other (i.e., bio2, mean diurnal temperature range; bio7, temperature annual range; bio13, precipitation of wettest month; bio15, precipitation seasonality). All these extra factors were previously interpolated to 100 m spatial resolution. Then, the final set of topoclimatic variables was composed by three topographic and four climatic layers.

Anthropic variables were defined taking into consideration the two sides of the spatial scale: grain (resolution) and extent (patch context). The variable “alteration” was built to quantify the direct habitat loss following anthropogenic alteration (sensu Fahrig [35]), and was modelled at the species scale by coding each cell with zero (not altered) or one (altered). Considering the historical background of land use in Sardinia, habitats completely not influenced by human activities are very scarce. Thus we defined “alteration” an intensive and direct exploitation of land, currently on-going and such that it prevents semi-natural structure from evolving (from this perspective, cropland and urban settlements are examples of altered habitat; cork-oak forest and pastures have to be considered not altered). We reclassified the land-use categories on the basis of this definition and the peculiar characteristic of Sardinian territory (table S2). To quantify the alteration of the context of each pixel at different spatial scales, we first chose a set of increasing distances, starting from the one that allow defining as “context” the eight pixels directly in contact with the considered one (150 m, 500 m, 1 km, 5 km, 10 km, 15 km, 20 km). For each distance we obtained the proportion of altered pixels (defined by variable “alteration”) within a circle centred in the centre of the pixel and with a radius equal to the distance considered. Then we computed the spatial correlation among all these “context alteration” layers (table S2) and we selected: i) the minimum distance at which the correlation coefficient (rs) with “alteration” did not exceed 0.70; ii) the subsequent minimum distances at which correlations with the previous selected layer fell below 0.70. In this way two distance ranges were chosen: 1 km (the finest scale for pixel context – variable PC-fine) and 15 km (variable PC-coarse). In conclusion, we derived three anthropic variables: “alteration” which represents the pixel state with respect to human disturbance; PC-fine and PC-coarse which measures context alteration at two different spatial scales. The spatial correlations among the final set of variables did not show coefficients exceeding 0.70 (maximum value was 0.68; table S3). Thus, we were able to build SDM under four causal hypotheses: H0 (topoclimatic), H1 (topoclimatic and alteration), H2 (topoclimatic, alteration and PC-fine) and H3 (topoclimatic, alteration and PC-coarse).

Modelling procedures

Many algorithms are available to model the spatial distribution of a species [6]. Among them, we chose the maximum entropy method [51] implemented by the software Maxent (version 3.3.3e; http://www.cs.princeton.edu/~schapire/maxent/), because it has been demonstrated to outperform other methods [52], and does not require absence data (it uses presence sites and background sites) [51]. When using Maxent, two fundamental decisions have to be taken. The first deals with the choice of the background. Besides the set of the presence points, indeed, Maxent uses a set of background points, which allow characterizing the environment available to the species. The choice of the background is particularly critical because it may influence the results if the sampling effort across the study area is not homogeneous [53]. A general framework to deal with sampling bias is getting the background points with the same bias as the presence sample [54]. Our sampled areas were not homogeneously distributed across the Sardinian island, but they concentrated along the road network (Figure S1). So, we decided to take the background points from a 2 km wide buffer along the roads (the map of the road network was obtained from http://www.sardegnageoportale.it). The selected distance represented the mean plus one standard deviation of the distances between each presence point and the nearest road. The background points were then generated by extracting 10,000 random points from this reduced area and adding the presence points. In this way almost the same sample bias affects presence sites and background points [54]. The second decision concerns the kind of “features” which have to be included in the model. Maxent, in fact, does not use the raw variables, but their transformed versions (called “features”) [51]. Maxent handles five kinds of features: linear, quadratic, product, threshold and hinge. The use of features instead of the raw variables is preferable when a complex process has to be modelled [55], but it leads to some over-fitting problems [55,56]. To avoid this theoretical risk we developed all the models by using only the hinge features, reducing redundancy [55].

Models evaluation and comparison

Model discriminating performance was estimated using AUC (Area Under the receiver operating characteristics Curve), which is a threshold-independent measure of the ability of the model to discriminate between background and presence sites [48]. AUC is widely applied in species distribution modelling to assess the classification ability of a model, but its employment in presence/background studies has been criticised by some authors [5759] because: i) “AUC assumes nothing about the relative costs of errors of omission and commission” [58]; ii) high AUC values may not reflect real accuracy if the “test presences disproportionally represent inhabited areas” [59]. In order to partially overcome these limits, we adopted two strategies: firstly, as suggested by Smith [59], we tested the significance of the observed AUC through randomization [60]. For both the species we randomly extracted from the background 999 samples of the same size as the observed ones (P. tiliguerta, N = 291; P. siculus, N = 355); we used these pseudoreplicates to generate 999 models and obtain the distribution of the AUC values. Then, we computed the probability associated to the observed AUC using the formula:

P=v+1n+1(eqn. 1)

where v is the number of values in the AUC distribution that exceed the observed one, and n is the total number of the pseudoreplicates [61]. Secondly, we converted the continuous suitability scores of Maxent output into binary outcomes (0 for absence, 1 for presence) and then we measured sensitivity (frequency of presence sites correctly classified) and specificity (frequency of pseudo-absences correctly classified) of each model. The cut-off value needed for the output conversion was chosen in order to maximize the two measures, i.e. to obtain equal sensitivity and specificity [62].

Model selection was carried out by means of IMTC [24]. We used the small sample unbiased Akaike Information Criterion (AICc) [63], since all the tested models had a number of free parameters that exceeded n/40 (n is the number of the observations; table 1) [32,33]. The models were ranked on the basis of the difference between the model AICc and the minimum value of AICc among the models. Also the Akaike weights (wi) were computed [33]. Akaike weights can be interpreted as the probability that model i is the best model for the observed data, given the candidate set of models [32]. The AICc were calculated with ENMtools ver. 1.3 (http://enmtools.blogspot.com/) [56,64].

SpeciesHyp.kAUCThres.Sens. = Spec.AICcΔi
P. siculusH2450.7740.3940.6926261.6800.000
H3460.7500.3910.6576314.61352.933
H0410.7450.3980.6656319.10557.425
H1440.7470.3950.6606319.51257.832
P. tiliguertaH2560.8110.4250.7315117.2080.000
H3540.7740.4530.7035194.40977.201
H1450.7680.4650.7005202.80785.599
>H0450.7470.4600.6785250.413133.205

Table 1. Models performance and comparison.

Hyp: hypothesis used to build the model (H0, topoclimatic; H1, topoclimatic + alteration; H2, topoclimatic + alteration + PC-fine; H3 topoclimatic + alteration + PC-coarse); k number of estimated parameters in the model. AUC: Area Under the ROC curve: probability values obtained by randomization tests are lower than 0.001 for all the models. Thres: threshold, cut-off value for binary conversion of the original continuous output. Sens. = Spec.: sensitivity and specificity; the two measures are the same because the threshold was chosen to maximize both of them. AICc: Akaike Information Criterion corrected for small sample size. Δi: AICc difference between model i and the best model. Akaike weights are not presented since they were one for the best model in both species.
CSV
Download CSV

Finally, in order to assess the effect size of the three hypothesized factors (topoclimatic, pixel and context alteration), we applied a variation partitioning analysis [65] adapted to univariate model output, following Muñoz et al. [66]: the approach allows disentangling the pure effect of each factor while controlling for between-factors overlap. To evaluate the direction of the effect, we visually investigated marginal and single response curves generated by Maxent: the former plots the change in model prediction as each variable is varied, keeping all the other ones at their average sample value; the latter is obtained by constructing a model using only one variable at a time, and then plotting the response with respect to the possible values of the variable itself [51]. The combination of the two ways allowed evaluating the effect of the among-variables correlations on the response. To obtain confidence intervals of the response curves, we used the cross-validation procedure of Maxent, splitting the original sample in ten sub-samples.

Results

A summary of the results concerning the developed SDM is reported in table 1. Looking at the discriminating performances (AUC), all the models show a significant deviation from chance and the best models for both species correspond to H2. Topoclimatic model (H0) is clearly the worst among those of P. tiliguerta, while it performs like H1 and H3 for P. siculus. The analysis of sensitivity and specificity values (table 1) leads to the same conclusion. The visual comparison of model maps (Figure 2) highlights that those built with PC-fine (H2) tend to better identify the unsuitable areas (blue colour), especially for the Tyrrhenian wall lizard. In both species, H2 hypothesis is also the most informative, showing the lowest AICc and having Akaike weight very close to one. The ranking of the models is slightly different between the two species (H2>H3>H1>H0 for P. tiliguerta; H2>H3>H0≥H1 for P. siculus), but all the models that have incorporated patch context characterizations are more informative than those that have simply used the alteration state of the pixels (H1).

thumbnail
Figure 2. Maps of the competing SDM.

A, B, C, D represent the models under hypotheses respectively H0, H1, H2 and H3 for P. siculus. E, F, G, H represent the same hypotheses for P. tiliguerta. Occurrences are also shown in separate maps.

https://doi.org/10.1371/journal.pone.0067573.g002

The variation partitioning of the output from H2 gives different results in the two lizards’ cases (Figure 3). The pure effect of climate and topography explains the largest part of model variation for P. siculus (80.61%), while the role of patch context is marginal (2.34%). Alteration shows nor a pure nor an overlapped effect on distribution and all the intersections are null. Results for P. tiliguerta are more complex: pure effects are found for climate and patch context (23.30% and 12.45% respectively) but not for alteration which shows complete overlap with patch context alone on one hand (6.32%), and with the other two factors together on the other hand (10.77%); topoclimatic and patch context factors share 11.67% of variation, after considering alteration.

thumbnail
Figure 3. Variation partitioning diagrams for the H2 hypothesis for the two species.

Circles represent variation explained by each factor (climate and topography, alteration, patch context). Numbers correspond to the percentage of variation associated to each circle subpart (pure, two-factor intersection, three factors intersection). The percentage associated to intersecting areas has not to be interpreted as interaction, but as a variation indifferently assignable to one or more factors [65]. Values smaller than 0.01% are not shown.

https://doi.org/10.1371/journal.pone.0067573.g003

The effects of the anthropogenic variables on the modelled suitability are graphically synthesized in Figure 4. All the response curves showed very narrow ranges, meaning low variation in response with the subsample considered. Marginal and single responses to habitat alteration show similar patterns: a negative relationship with suitability for the Tyrrhenian wall lizard, weak in the marginal response (Figure 4A), more evident in the single response (Figure 4B); no apparent relationship for the Italian wall lizard (Figure 4A, 4B). Also the responses to patch context (PC-fine) are quite different between the two species (Figure 4C, 4D): the Italian wall lizard seems to tolerate higher levels of altered context than the congeneric species, and suitability starts to decrease when alteration is above 40% (marginal response) or 85% (single response). On the contrary, P. tiliguerta appears to be more sensitive to altered context, since, in both graphs, suitability decreases with increasing alteration: the reduction becomes faster when alteration is above 40% (Figure 4C, 4D).

thumbnail
Figure 4. Response curve for the anthropogenic variable of the best model (H2).

A: marginal response curves for the variable “alteration” (suitability is obtained from the full model keeping constant all the variables but alteration). B: single variable response curves for the same variable (suitability is obtained from a model including only the variable “alteration”). C: marginal response curves for PC-fine (patch context alteration). D: single variable response curves for the same variable. Solid lines indicate the mean of ten cross validated models; dashed lines represent the minimum and maximum range of the response curves.

https://doi.org/10.1371/journal.pone.0067573.g004

Discussion

During recent years, scale dependent effects are attracting more and more interests of the studies on species distribution modelling [714]. There is increasing evidence that reliability of SDM as well as our understanding of the consequences of human-induced alterations largely depend on the scale at which they are modelled. This may have relevant effects on the conservation measures that can be implemented to preserve species. We proposed a generalizable two-steps approach to study the effect of the anthropogenic alteration of the patch context [19] on the spatial distribution of a species, in order to define up to what extent a distance can be regarded as “context” for a given species in a given geographical region. In our exemplification, based on the occurrence data of two lacertid lizards, the models built with the inclusion of variables accounting for habitat alteration of the context (PC-fine, PC-coarse) have proven to perform better than the basal models (table 1), and, more interestingly, than the models which incorporate only local alteration (H1). For both species, the best results were obtained setting the extent of patch context to 1 km (PC-fine, H2; table 1), and this distance has produced the strongest trimming effect on the potential distribution (Figure 2). PC-fine may be seen as a measure of isolation of each cell [35,67], a property strictly related to accessibility, the “M" component of the BAM diagram [17] [18],: the permeability of the context towards movements decreases as isolation increases, thus preventing to access the suitable patches [3] and negatively interfering with meta-population dynamics [68]. Clearly, only a variable that works at the context level of the spatial scale is able to collect information about patch isolation and consequently its use in model building is expected to produce more realistic SDM. Now the question becomes: how far do we have to go by the pixel? Our results show that an increase in the extent of the context from 1 to 15 km produces less informative model, thus indicating that a 15 km range is too large. The outcome seems plausible considering that M-factor does not depend only on the permeability of the context, but also on the dispersal ability of the species. In the case study this ability was about 100 m [4649], so the extent of the context is not expected to become too large before a reduction in its effect on shaping the species distribution may occur. Previous studies about the influence of patch context on animals with relative small dispersal ability (but which made use of data from well surveyed sites and of other inference techniques), have set the effect at the scale of some kilometres (reptile and amphibian: 2 km [69]; beetles: till 2 km [19]; butterflies: from 1 km to 4 km [9]).

An additional relevant issue of this study is that the alteration of the pixel context has larger consequences than the alteration of the pixel itself (Figure 3). Also this result can be interpreted in the light of the study species and the definition of anthropogenic alteration we employed. Actually, both lizards (even though differences occur) tolerate human proximity, and may inhabit some kinds of anthropogenic habitats [40,43]. In this case suitability depends on the “semi-natural” microhabitats (e.g., gardens, rows of trees, ecotones between crops, wooded micro-patches) available inside each 100 × 100 m cell. The degree of suitability guaranteed by the occurrence of these micro-habitats depends on their amount and quality in addition to the species considered (e.g., the endemic lizard is more exacting and sensitive to intensive land exploitation [39]). Nevertheless, the association of contiguous altered patches may become a limiting factor for dispersal when there is a lack of continuity in suitable microhabitats among adjacent cells: this condition is able to shape the species distribution through the M-factor. The above interpretation is consistent with the analysis of the variation partitioning results (Figure 3) and the response curves (Figure 4), which show that: i) alteration of the pixel does not substantially affect suitability (no pure effect were found in both species) and the apparent trend exhibited in the single response graph of P. tiliguerta is probably due to a correlation with PC-fine; ii) altered context (PC-fine) heavily penalizes P. tiliguerta, whereas it is partly tolerated by P. siculus. This is what we might expect considering the different characteristics of the two species. In fact, the Italian wall lizard: (i) is more synanthropic than the congeneric Tyrrhenian wall lizard [39,42]; (ii) is more effective in thermoregulation (its body temperature shows little variation when the environmental temperature changes [42]), so it stands up well to the warm conditions typical of the anthropic habitats; (iii) is a better performer in locomotor endurance tests [70]; (iv) is highly effective in colonizing anthropic habitats, as long as parks and garden are available [38]. All these characteristics lead us to claim that P. siculus should have better dispersal ability in altered landscapes, and that isolation should be a less effective obstacle for this species, at least up to a certain threshold.

From the conservation point of view, our case study allows drawing some important guidelines, which may be the starting point for management actions and/or for further investigation. Firstly, the endemic species would be seriously threatened by habitat alteration rising; the same would not be true for the introduced as well as long time naturalized species. So if other natural areas will be lost, we might assist to a regression of the endemic lizard, without a meaningful change in the distribution of the introduced P. siculus. Secondly, the main cause of concern at current state is represented by habitat isolation, measured at one kilometre scale; so it would be worthwhile to focus any actions primarily to this scale in order to limit the rise of this parameter value. The importance of establishing the scale and the aim of a conservation plan does not need further discussion [15,16].

Generalizing our findings, the proposed approach, which aimed to link the spatial distribution of a species to the anthropogenic alteration of the patch context, has proved to work well, giving reasonable and ecologically interpretable results. It combines the advantages of SDM and of ITMC. The former allows using available data on species distribution, not necessarily taken for conservation aims [71], so reducing the field work effort and the related costs. Furthermore it allows simultaneously analysing larger portion of the distribution area of a species than traditional approaches, which typically use few subsamples of it [9,13,15,23]. The latter (ITMC), being a general method to compare explicit hypotheses and weigh quantifiable variable importance [32,33], may be easily applied to assess the scale at which human disturbance produces the strongest effects. The same approach can work well even with other measurement of human disturbance (for instance: pollutant concentration, invasive species presence, pesticides, road network extension): it is sufficient to translate human disturbance into variables quantifiable at different scales and then compare the models built under the null and the alternative hypotheses.

Finally, we underline that a critical step in applying our procedure is the choice of the “right” investigating scale (i.e., the resolution of the environmental and occurrence data used in developing the SDM). This choice is fundamental in order to correctly interpret patch context effect: using the same order of magnitude of the species dispersal ability may be a good idea [67], because it allows disentangling the effects of habitat loss from those of fragmentation/isolation at the species point of view [15].

Supporting Information

Figure S1. Map of the occurrence points and the road network.

https://doi.org/10.1371/journal.pone.0067573.s001

(DOC)

Table S1. Spatial Correlations among topographic and climatic variables.

https://doi.org/10.1371/journal.pone.0067573.s002

(DOC)

Table S3. Spatial correlations among the final candidate variables.

https://doi.org/10.1371/journal.pone.0067573.s004

(DOC)

Acknowledgments

All the colleagues and in particular Carmen Fresi and Maria Grazia Satta involved in fieldwork. We also thank Diego Fontaneto and four anonymous reviewers for their valuable comments that sensibly improved an early draft of the article.

Author Contributions

Conceived and designed the experiments: MM RS SS. Performed the experiments: LB VN CC. Analyzed the data: MM RS SS. Contributed reagents/materials/analysis tools: MM RS SS. Wrote the manuscript: MM.

References

  1. 1. Sala OE, Chapin FS, Armesto JJ, Berlow E, Bloomfield J et al. (2000) Global biodiversity scenarios for the year 2100. Science 287: 1770-1774. doi:https://doi.org/10.1126/science.287.5459.1770. PubMed: 10710299.
  2. 2. Wilson EO (2002) The future of life. London: Little. Brown.: 256.
  3. 3. Peterson AT, Soberón J (2012) Species distribution modeling and ecological niche modeling: getting the concepts right. Nat Conservacao 10: 1-6. doi:https://doi.org/10.4322/natcon.2012.001.
  4. 4. Guisan A, Thuiller W (2005) Predicting species distribution: offering more than simple habitat models. Ecol Lett 8: 993-1009. doi:https://doi.org/10.1111/j.1461-0248.2005.00792.x.
  5. 5. Pearson RG (2010) Species’ Distribution Modeling for Conservation Educators and Practitioners. Lessons Conserv 3: 54-89.
  6. 6. Peterson AT, Soberón J, Pearson RG, Anderson RP, Martínez-Meyer E et al. (2011) Ecological niches and geographic distributions. Princeton: Princeton University Press. 328pp.
  7. 7. Graf RF, Bollmann K, Suter W, Bugmann H (2005) The importance of spatial scale in habitat models: capercaillie in the Swiss Alps. Landscape Ecol 20: 703-717. doi:https://doi.org/10.1007/s10980-005-0063-7.
  8. 8. Guisan A, Graham CH, Elith J, Huettmann F, Dudik M et al. (2007) Sensitivity of predictive species distribution models to change in grain size. Divers Distrib 13: 332-340. doi:https://doi.org/10.1111/j.1472-4642.2007.00342.x.
  9. 9. Cozzi G, Muller CB, Krauss J (2008) How do local habitat management and landscape structure at different spatial scales affect fritillary butterfly distribution on fragmented wetlands? Landscape Ecol 23: 269-283. doi:https://doi.org/10.1007/s10980-007-9178-3.
  10. 10. Menke SB, Holway DA, Fisher RN, Jetz W (2009) Characterizing and predicting species distributions across environments and scales: Argentine ant occurrences in the eye of the beholder. Glob Ecol Biogeogr 18: 50-63. doi:https://doi.org/10.1111/j.1466-8238.2008.00420.x.
  11. 11. Hortal J, Roura-Pascual N, Sanders NJ, Rahbek C (2010) Understanding (insect) species distributions across spatial scales. Ecography 33: 51-53. doi:https://doi.org/10.1111/j.1600-0587.2009.06428.x.
  12. 12. Convertino M, Kiker GA, Muñoz-Carpena R, Chu-Agor ML, Fischer Ra et al . (2011) Scale- and resolution-invariance of suitable geographic range for shorebird metapopulations. Ecol Complex 8: 364-376. doi:https://doi.org/10.1016/j.ecocom.2011.07.007.
  13. 13. Martin AE, Fahrig L (2012) Measuring and selecting scales of effect for landscape predictors in species-habitat models. Ecol Appl 22: 2277-2292. doi:https://doi.org/10.1890/11-2224.1. PubMed: 23387125.
  14. 14. Pineda E, Lobo JM (2012) The performance of range maps and species distribution models representing the geographic variation of species richness at different resolutions. Glob Ecol Biogeogr 21: 935-944. doi:https://doi.org/10.1111/j.1466-8238.2011.00741.x.
  15. 15. Wiens JA (1989) Spatial Scaling in Ecology. Funct Ecol 3: 385-397. doi:https://doi.org/10.2307/2389612.
  16. 16. Levine SA (1992) The Problem of Pattern and Scale in Ecology: The Robert H MacArthur Award Lecture. Ecology 73: 1943-1967. doi:https://doi.org/10.2307/1941447.
  17. 17. Saupe EE, Barve V, Myers CE, Soberón J, Barve N et al. (2012) Variation in niche and distribution model performance: the need for a priori assessment of key causal factors. Ecol Modell 237-238: 11-22. doi:https://doi.org/10.1016/j.ecolmodel.2012.04.001.
  18. 18. McGill BJ (2010) Ecology. Matters of scale. Science 328: 575–576. doi:https://doi.org/10.1126/science.1188528. PubMed: 20431001.
  19. 19. Holland JD, Bert DG, Fahrig L (2004) Determining the spatial scale of species’ response to habitat. BioScience 54: 227-233.
  20. 20. Wiens JA, Stenseth NC, Van Horne B, Ims RA (1993) Ecological mechanisms and landscape ecology. Oikos 66: 369-380. doi:https://doi.org/10.2307/3544931.
  21. 21. Mazerolle MJ, Villard M (1999) Patch characteristics and landscape context as predictors of species presence and abundance: a review. EcoScience 6: 117-124.
  22. 22. Debinski DM, Ray C, Saveraid EH (2001) Species diversity and the scale of the landscape mosaic: do scales of movement and patch size affect diversity? Biol Conserv 98: 179-190. doi:https://doi.org/10.1016/S0006-3207(00)00153-1.
  23. 23. Kattwinkel M, Strauss B, Biedermann R, Kleyer M (2009) Modelling multi-species response to landscape dynamics: mosaic cycles support urban biodiversity. Landscape Ecol 24: 929-941. doi:https://doi.org/10.1007/s10980-009-9371-7.
  24. 24. Stephens PS, Buskirk SW, Hayward GD, Martinez del Rio C (2005) Information theory and hypothesis testing: a call for pluralism. J Appl Ecol 42: 4-12. doi:https://doi.org/10.1111/j.1365-2664.2005.01002.x.
  25. 25. McInerny GJ, Etienne RS (2012) Ditch the niche – is the niche a useful concept in ecology or species distribution modelling? J Biogeogr 39: 2096-2102. doi:https://doi.org/10.1111/jbi.12033.
  26. 26. McInerny GJ, Etienne RS (2012) Stitch the niche –a practical philosophy and visual schematic for the niche concept. J Biogeogr 39: 2103-2111. doi:https://doi.org/10.1111/jbi.12032.
  27. 27. McInerny GJ, Etienne RS (2012) Pitch the niche – taking responsibility for the concepts we use in ecology and species distribution modelling. J Biogeogr 39: 2112-2118. doi:https://doi.org/10.1111/jbi.12031.
  28. 28. Elith J, Leathwick JR (2009) Species Distribution Models: Ecological Explanation and Prediction Across Space and Time. Annu. Rev Ecol Evol S 40: 677-697. doi:https://doi.org/10.1146/annurev.ecolsys.110308.120159.
  29. 29. Brito JC, Fahd S, Martinez-Freiria F, Tarroso P, Larbes S et al. (2011) Climate change and peripheral populations: predictions for a relict Mediterranean viper. Acta Herpetol 6: 105-118.
  30. 30. Sillero N (2009) Modelling a species in expansion at local scale: is Hyla meridionalis colonising new areas in Salamanca, Spain? Acta Herpetol 4: 37-46.
  31. 31. Scali S, Mangiacotti M, Sacchi R, Gentilli A (2011) A tribute to Hubert Saint Girons: niche separation between Vipera aspis and V. berus on the basis of distribution models. Amphibia Reptilia 32: 223-233. doi:https://doi.org/10.1163/017353711X562171.
  32. 32. Johnson JB, Omland KS (2004) Model selection in ecology and evolution. Trends Ecol Evol 19: 101-108. doi:https://doi.org/10.1016/j.tree.2003.10.013. PubMed: 16701236.
  33. 33. Anderson DR, Burnham KP, Thompson WL (2000) Null Hypothesis Testing: Problems, Prevalence, and an Alternative. J Wildl Manage 64: 912-923. doi:https://doi.org/10.2307/3803199.
  34. 34. Acevedo P, Jiménez-Valverde A, Lobo JM, Real R (2012) Delimiting the geographical background in species distribution modelling. J Biogeogr 39: 1383-1390. doi:https://doi.org/10.1111/j.1365-2699.2012.02713.x.
  35. 35. Fahrig L (2003) Effects of habitat fragmentation on biodiversity. Annu. Rev Ecol Evol S 34: 487-515. doi:https://doi.org/10.1146/annurev.ecolsys.34.011802.132419.
  36. 36. Peel MC, Finlayson BL, McMahon TA (2007) Updated world map of the Köppen-Geiger climate classification. Hydrol Earth Syst Sci 11: 1633-1644. doi:https://doi.org/10.5194/hess-11-1633-2007.
  37. 37. Bacchetta G, Bagella S, Biondi E, Farris E, Filigheddu RS et al. (2009) Vegetazione forestale e serie di vegetazione della Sardegna (con rappresentazione cartografica alla scala 1:350.000). Fitosociologia 46, Suppl. 1: 3-82.
  38. 38. Sindaco R, Doria G, Razzetti E, Bernini F (2006) Atlante degli anfibi e dei rettili d’Italia / Atlas of Italian amphibians and reptiles. Firenze: Societas Herpetologica Italica. Edizioni Polistampa. 800pp.
  39. 39. Corti C, Lo Cascio P (2002) The lizards of Italy and adjacent areas. Frankfurt am Main: Edition Chimaira. 165 pp.
  40. 40. Bruschi S, Corti C, Capula M (2010) Podarcis tiliguerta (Gmelin, 1789). In: C. CortiM. CapulaL. LuiselliE. RazzettiR. Sindaco. Fauna d’Italia. Reptilia. Bologna: Edizioni Calderini de Edizioni Calderini de Il Sole 24 Ore Editoria Specializzata S.r.l.. pp 417-427.
  41. 41. Delaugerre M, Cheylan M (1992) Atlas de repartition des Batraciens et Reptiles de Corse. Pampelune: Parc natural regional de la Corse. École Pratique des Hautes Etudes. 128pp.
  42. 42. Van Damme R, Bauwens D, Castilla AM, Verheyen R (1990) Comparative thermal ecology of the sympatric lizards Podarcis tiliguerta and Podarcis sicula. Acta Oecol 11: 503-512.
  43. 43. Corti C, Biaggini M, Capula M (2010) Podarcis siculus (Rafinesque-Schmaltz, 1810). In: C. CortiM. CapulaL. LuiselliE. RazzettiR. Sindaco. Fauna d’Italia. Reptilia 24 24 24. Bologna: Edizioni Calderini de Il Sole Ore Editoria Specializzata S.r.l. pp 407-417.
  44. 44. Corti C, Böhme W, Delfino M, Massetti M (1999) Man and lacertids on the Mediterranean islands: conservation perspectives. Nat Croat 8: 287-300.
  45. 45. Podnar M, Mayer W, Tvrtković N (2005) Phylogeography of the Italian wall lizard, Podarcis sicula, as revealed by mitochondrial DNA sequences. Mol Ecol 14: 575-588. doi:https://doi.org/10.1111/j.1365-294X.2005.02427.x. PubMed: 15660947.
  46. 46. Foà A, Bearzi M, Baldaccini NE (1990) Preliminary report on the size of the home range and on the orientational capabilities in the lacertid lizard Podarcis sicula. Ethol Ecol Evol 2: 310. doi:https://doi.org/10.1080/08927014.1990.9525443.
  47. 47. Brown RM, Gist DH, Taylor DH (1995) Home range ecology of an introduced population of the European wall lizard Podarcis muralis (Lacertilia; Lacertidae) in Cincinnati, Ohio. Am Midl Nat 133: 344-359. doi:https://doi.org/10.2307/2426399.
  48. 48. Perry G, Garland T Jr (2002) Lizard home ranges revisited: effects of sex, body size, diet, habitat, and phylogeny. Ecology 83: 1870-1885. doi:https://doi.org/10.1890/0012-9658(2002)083[1870:LHRREO]2.0.CO;2.
  49. 49. Scali S, Sacchi R, Azzusi M, Daverio S, Oppedisano T et al. (2013) Homeward bound: factors affecting homing ability in a polymorphic lizard. J Zool 289: 296-203.
  50. 50. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A (2005) Very high resolution interpolated cli-mate surfaces for global land areas. Int J Climatol 25: 1965-1978. doi:https://doi.org/10.1002/joc.1276.
  51. 51. Phillips SJ, Anderson RP, Schapire RE (2006) Maximum entropy modeling of species geographic distributions. Ecol Modell 190: 231-259. doi:https://doi.org/10.1016/j.ecolmodel.2005.03.026.
  52. 52. Elith J, Graham CH, Anderson RP, Dudík M, Ferrier S et al. (2006) Novel methods improve prediction of species’ distributions from occurrence data. Ecography 29: 129-151. doi:https://doi.org/10.1111/j.2006.0906-7590.04596.x.
  53. 53. Barbet-Massim M, Jiguet F, Albert CH, Thuiller W (2012) Selecting pseudo-absences for species distribution models: how, where and how many? Methods Ecol Evol 3: 327-338. doi:https://doi.org/10.1111/j.2041-210X.2011.00172.x.
  54. 54. Phillips SJ, Dudík M, Elith J, Graham CH, Lehmann A et al. (2009) Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data. Ecol Appl 19: 181-197. doi:https://doi.org/10.1890/07-2153.1. PubMed: 19323182.
  55. 55. Elith J, Phillips SJ, Hastie T, Dudík M, Chee YE et al. (2010) A statistical explanation of MaxEnt for ecologists. Divers Distrib 17: 43-57.
  56. 56. Warren DL, Seifert SN (2011) Ecological niche modeling in Maxent: the importance of model complexity and the performance of model selection criteria. Ecol Appl 21: 335-342. doi:https://doi.org/10.1890/10-1171.1. PubMed: 21563566.
  57. 57. Peterson AT, Papeş M, Soberón J (2008) Rethinking receiver operating characteristic analysis applications in ecological niche modeling. Ecol Modell 213: 63-72. doi:https://doi.org/10.1016/j.ecolmodel.2007.11.008.
  58. 58. Jiménez-Valverde A (2012) Insights into the area under the receiver operating characteristic curve (AUC) as a discrimination measure in species distribution modelling. Glob Ecol Biogeogr 21: 498-507. doi:https://doi.org/10.1111/j.1466-8238.2011.00683.x.
  59. 59. Smith AB (2013) On evaluating species distribution models with random background sites in place of absences when test presences disproportionately sample suitable habitat. Divers Distrib. doi:https://doi.org/10.1111/ddi.12031.
  60. 60. Raes N, ter Steege H (2007) A null-model for significance testing of presence-only species distribution models. Ecography 30: 727-736. doi:https://doi.org/10.1111/j.2007.0906-7590.05041.x.
  61. 61. Hestelberg T, Monaghan S, Moore DS, Clipson A, Epstein R (2003) Bootstrap methods and permutation tests. New York: W. H. Freeman and Company. 81pp.
  62. 62. Freeman EA, Moisen GG (2008) A comparison of the performance of threshold criteria for binary classification in terms of predicted prevalence and kappa. Ecol Modell 217: 48-58. doi:https://doi.org/10.1016/j.ecolmodel.2008.05.015.
  63. 63. Akaike H (1973) Information theory as an extension of the maximum likelihood principle. In: BN PetrovF. Csaki. Second International Symposium on Information Theory. . Budapest: Akadémiai Kiadó. pp. 267-281.
  64. 64. Warren DL, Glor RE, Turelli M (2010) ENMTools: a toolbox for comparative studies of environmental niche models. Ecography 33: 607-611.
  65. 65. Legendre P, Legendre L (1998) Numerical ecology. Second English edition. Amsterdam: Elsevier Science B.V.. p. 870.
  66. 66. Muñoz AR, Real R, Barbosa AM, Vargas JM (2005) Modelling the distribution of Bonelli’s eagle in Spain: implications for conservation planning. Divers Distrib 11: 477-486. doi:https://doi.org/10.1111/j.1366-9516.2005.00188.x.
  67. 67. Bender DJ, Tischendorf L, Fahrig L (2003) Using patch isolation metrics to predict animal movement in binary landscapes. Landscape Ecol 18: 17-39. doi:https://doi.org/10.1023/A:1022937226820.
  68. 68. Hanski I (1998) Metapopulation dynamics. Nature 396: 41-49. doi:https://doi.org/10.1038/23876.
  69. 69. Findlay CS, Houlahan J (1997) Anthropogenic correlates of species richness in southeastern Ontario wetlands. Conserv Biol 11: 1000-1009. doi:https://doi.org/10.1046/j.1523-1739.1997.96144.x.
  70. 70. Vanhooydonk B, Van Damme R, Aerts P (2000) Ecomorphological correlates of habitat partitioning in Corsican lacertid lizards. Funct Ecol 14: 358-368. doi:https://doi.org/10.1046/j.1365-2435.2000.00430.x.
  71. 71. Smith AB (2013) The relative influence of temperature, moisture and their interaction on range limits of mammals over the past century. Glob Ecol Biogeogr 22: 334-343. doi:https://doi.org/10.1111/j.1466-8238.2012.00785.x.