Incorporating natural and human factors in habitat modelling and spatial prioritisation for the Lynx lynx martinoi

Countries in south-eastern Europe are cooperating to conserve a sub-endemic lynx species, Lynx lynx martinoi. Yet, the planning of species conservation should go hand-in-hand with the planning and management of (new) protected areas. Lynx lynx martinoi has a small, fragmented distribution with a small total population size and an endangered population. This study combines species distribution modelling with spatial prioritisation techniques to identify conservation areas for Lynx lynx martinoi. The aim was to determine locations of high probability of occurrence for the lynx, to potentially increase current protected areas by 20 % in Albania, the former Yugoslav Republic of Macedonia, Montenegro, and Kosovo. The species distribution modelling used generalised linear models with lynx occurrence and pseudo-absence data. Two models were developed and fitted using the lynx data: one based on natural factors, and the second based on factors associated with human disturbance. The Zonation conservation planning software was then used to undertake spatial prioritisations of the landscape using the first model composed of natural factors as a biological feature, and (inverted) a second model composed of anthropological factors such as a cost layer. The first model included environmental factors as elevation, terrain ruggedness index, woodland and shrub land, and food factor as chamois prey (occurrences) and had a prediction accuracy of 82 %. Second model included anthropological factors as agricultural land and had a prediction accuracy of 65 %. Prioritised areas for extending protected areas for lynx conservation were found primarily in the Albania–Macedonia–Kosovo and Montenegro–Albania–Kosovo cross-border areas. We show how natural and human factors can be incorporated into spatially prioritising conservation areas on a landscape level. Our results show the importance of expanding the existing protected areas in cross-border areas of core lynx habitat. The priority of these cross-border areas highlight the importance international cooperation can play in designing and implementing a coherent and long-term conservation plan including a species conservation plan to securing the future of the lynx. Published by Copernicus Publications on behalf of the European Ecological Federation (EEF). 18 K. Laze and A. Gordon: Spatial prioritisation for Lynx lynx martinoi


Introduction
Many large carnivore species have declining population levels, and are listed as endangered in the countries they inhabit (Wiegand et al., 1998;Fernández et al., 2003;Kramer-Schadt et al., 2005). A prime example is Lynx lynx martinoi that is one of nine subspecies (namely Lynx, carpathicus, martinoi, dinniki, isabellinus, wardi, kozlovi, wrangeli, stroganovi) distributed in Europe, central Asia, and Russia (Arx et al., 2004). In 2005, Lynx lynx martinoi was shown to be a sub-species of the Eurasian lynx (Lynx lynx) using genetic analysis (Breitenmoser-Würsten and Breitenmoser, 2001). The Lynx lynx martinoi once populated a large area stretching from Slovenia, Serbia and Croatia to Albania, the former Yugoslav Republic of Macedonia, Montenegro, and Kosovo (Bojovič, 1978;Mirič, 1978; Large Carnivore Initiative for Europe, 1997Europe, , 2004Breitenmoser et al., 2008). Today, the Balkan Lynx Recovery Program is contributing by collecting records of the occurrences of individuals and populations of Lynx lynx martinoi in Albania, the former Yugoslav Republic of Macedonia, Kosovo, and Montenegro (IUCN/SSC Cat Specialist Group, 2011;Balkan Lynx Recovery Programme, 2015). The Balkan Lynx Recovery Programme has done work analysing the distribution of the lynx population, the human factors impacting its distribution and drawing and implementing strategies on the conservation of the lynx (Balkan Lynx Strategy Group, 2008). Explanatory factors impacting the spatial distribution of the Lynx lynx martinoi, and the prioritisation of areas in the landscape for the conservation are less known and studied than other subspecies. The latest population estimates is below 50 mature individuals, resulting in the Lynx lynx martinoi being listed as critically endangered on the IUCN Red List of Threatened Species on 19 November 2015 (IUCN, 2015).
The Lynx lynx martinoi (henceforth referred to as "the lynx") occupies deciduous forests containing tree species of Fagus sylvatica, Quercus spp., Carpinus betulus, Ostrya carpinifolia, Fraxinus ornus. The lynx occupies coniferous forests of Pinus spp., of Abies spp., as well as mixed forests of fir and beech. It prefers elevated areas and uses mountain pastures for hunting in summer. Preliminary results obtained from a study of radio-telemetry in the Former Yugoslav Republic (FYR) of Macedonia revealed that the lynx diet comprised of chamois (Rupicapra rupicapra; 24 %), European brown hare (Lepus europaeus; 12 %) and roe deer (Capreolus capreolus; 64 %) (IUCN/SSC Cat Specialist Group, 2012).
Today, the most important threats to the lynx are a small population (less than 50 individuals, increasing the risk of extinction from demographic stochasticity), a limited prey base (caused by illegal hunting in Albania), degradation of habitat (particularly in Albania and Kosovo), poaching activities by humans and the isolation of smaller populations in a fragmented landscape (IUCN/SSC Cat Specialist Group, 2012).
The viability of a species in the landscape is fundamentally linked with identification and delineation of its habitat. This is important for determining both the distribution and abundance of the species as well as for prioritising conservation activities targeting the species (Boyce et al., 2016). Species distribution models (SDMs) are widely used to identify where habitat for a species is likely to occur and to determine the core areas important for the conservation of species (Zielinski et al., 2006). The advantage of SDMs is that based on a given number of species records, the model can then make predictions on large scales as to where the habitat for a species is likely to occur. SDMs are increasingly using larger and more complex data sets of species due to greater amounts of data being collected (e.g. radio telemetry), the increased availability of environmental data from remote sensing and of powerful techniques like generalised linear models and geographical information system (GIS) to quantify species-environment relationships . Today, SDMs are the main tools used to produce spatially explicit predictions regarding the relationships between a species and its environment (Elith and Leathwick, 2009) at different spatial scales ranging from the local scale (e.g.  to the global scale (e.g. Thuiller et al., 2005). SDMs for one or more species are now commonly used as inputs for spatial conservation prioritisation, often focused on determining where new conservation areas should occur in the landscape (Moilanen et al., 2006) using the software of Zonation (Moilanen et al., 2012). Zonation software is an approach to reduce uncertainties of SDM results (Guisan et al., 2013) that may potentially be used by nature conservation decision making. SDM applications to support nature conservation decision making is still needed, although critical habitat and reserve selection are two of conservation domains where SDMs can be increasingly valuable (see Guisan et al., 2013).
Species present in existing protected areas are more likely to have viable populations if the areas surrounding them are protected (see, e.g. Fischer et al., 2006). Today, the conservation efforts of the lynx are focused on collecting data for existing population of the species, and on increasing the awareness of local human populations  of the existence of the species (Balkan Lynx Strategy Group, 2008). In addition, efforts have focused on gathering data on species reproduction occurrences and on identifying habitat for the species (e.g. in Montenegro, Balkan Lynx Recovery Programme, 2015). Yet to our knowledge, there are no studies on the following: habitat selection using a resource selection function (Boyce et al., 2016), developing a SDM, on animal behaviour and animal movement (Patterson et al., 2008), or on conservation prioritisation (Moilanen et al., 2006) for the lynx. A species distribution modelling and conservation prioritisation approach can help first identify the entire area of predicted habitat for the lynx to support future data collection for the species and to identify core conservation areas for the species that should be protected and priority areas for managing given existing protected areas. We provided the first large-scale SDM for the lynx to estimate the prob-ability of occurrence of the species throughout its range. The study area covered the full range of the current habitat of the lynx in south-eastern Europe. We applied species distribution modelling to understand how the lynx used natural resources (forests, prey). Species distribution modelling was combined with a spatial prioritisation approach to address questions such as the most effective places to extend existing conservation areas  for the lynx.
We aimed to develop SDM as a "model 1" composed of natural factors and "model 2" composed of anthropological factors: (1) to identify the areas of estimated probability of occurrence for the lynx and natural factors that explained the occurrence (records) of the lynx; (2) to understand the impact of anthropogenic land-use activities on the lynx distribution. The third aim was to identify priority areas for lynx conservation and to determine the locations for the expansion of existing protected areas in the study area.
Finally, we used the lynx model 1 (selected from species distribution modelling) derived from natural factors as a biological feature layer and the lynx model 2 (selected from species distribution modelling) derived from anthropogenic land-use activities as cost layer, respectively, in Zonation conservation planning software (Moilanen et al., 2012).

Occurrence records of lynx and lynx prey
Lynx data were collected by two domestic non-governmental organisations (Sect. S0 in the Supplement) using a 50question face-to-face interview with key informants (i.e. hunters, game wardens, foresters, shepherds, livestock breeders, beekeepers, cafeteria or market owners) in randomly selected villages in the northern and eastern Albania (91 villages) and western the FYR of Macedonia (154 villages) (Sect. S0) . At least two people were randomly selected to be interviewed in each village by Ivanov et al. (2008). A grid cell of 10 km × 10 km for Albania and the FYR of Macedonia (henceforth "Macedonia") was used based on the lynx population and lynx population density in 2001 of the publication of Arx et al. (2004)  . The population density of the lynx was estimated to be between 0.65 and 1.09 in Albania and 2.06 in Macedonia in 2001 (see Sect. S0). The spatial coverage of this survey was approximately 13 600 km 2 (i.e. 63 grid cells in Albania and 73 grid cells in Macedonia, Ivanov et al., 2008). The occurrence of the lynx and lynx prey (chamois, brown hare, roe deer) was based on the "relative number of positive answers" of interviewees on the species occurrences (Sect. 1 of the questionnaire was on presence and distribution of large mammal species in the last 5 years; Sect. S0), in a 10 km × 10 km grid cell. A probable occurrence (permanent occurrences in this study) was selected if there were more than 50 % affirmative answers (to questions, e.g. on hard facts like lynx tracks, stuffed lynx, prepared lynx pelts, lynx attacks on humans, lynx attacks on domestic animals, more than one lynx or female with cubs observed) of interviewees on the species occurrence in a grid cell of 10 km × 10 km . A possible occurrence (temporal occurrences in this study) resulted if there were less than 50 % of affirmative answers of interviewees on the species occurrence in a grid cell of 10 km × 10 km and no species occurrence if there were no affirmative answers (i.e. on hard facts like lynx tracks, stuffed lynx, prepared lynx pelts, lynx attacks on humans, lynx attacks on domestic animals, more than one lynx or female with cubs observed) of the interviewees on the species occurrences in a grid cell of 10 km × 10 km . See Sect. S0, for detailed information on lynx data collected between 2006 and 2007, as well as for data on lynx distribution and lynx data in 2001.
The Balkan Lynx Recovery Programme collected the lynx prey occurrence records consisting of species of brown hare, chamois and roe deer . Roe deer occurrences were located mostly in the same grid cell of 10 km × 10 km of the lynx occurrences . Brown hare and chamois had the widest and the smallest distribution, respectively. Chamois and brown hare records were used in this study to investigate any spatial effects of the smallest distribution of chamois and the widest distribution of brown hare in the occurrence (records) of the lynx in terms of neighbourhood scale of the best-performing model 1 (see Model selection for the lynx). In total, 87 % of records of chamois and 89 % of records of brown hare were located from approximately 2 to 23 km from the closest lynx occurrences. For more information on lynx prey occurrence data collected between 2006 and 2007, see Sect. S0.
We used 109 (39 permanent and 70 temporal) lynx occurrence records, 114 chamois occurrence records, 135 occurrence brown hare records in protected areas and public-owned land, in our study area. The lynx occurrence records (109) consisted of 37 lynx permanent occurrences records (25 in Macedonia and 12 in Albania) and 72 lynx temporal occurrences records (36 in Macedonia, 26 in Albania, 6 in Montenegro, and 4 in Kosovo). We geo-referenced lynx occurrence records from for the lynx in Kosovo and Montenegro from the maps of Arx et al. (2004), and lynx occurrence records from 2006 to 2007 in Albania and Macedonia from the maps of Ivanov et al. (2008) (Fig. 1). We also geo-referenced occurrence records of species of lynx prey of a resolution 10 km × 10 km for Albania and Macedonia from the maps of Ivanov et al. (2008), the Balkan Lynx compendium (IUCN/SSC Cat Specialist Group, 2011). We then obtained the coordinates of X and Y of the lynx occurrences and of the species of lynx prey (chamois and brown hare). We assigned a location for each species (lynx, chamois, brown hare) within a 1 km × 1 km grid cell study unit by spatially joining the layer of the species occurrence record and of the grid 1 km × 1 km using Spatial Join of the Analysis Tool in ArcGIS transferring the attribute of species occurrence record (X, Y ) to the 1 km × 1 km grid (for each species i.e. lynx, chamois, brown hare).
We argue that the lynx may be anywhere in a 10 × 10 km cell; however, the lynx is also thought to perceive forests up to 1 km apart as connected (Kramer-Schadt et al., 2005). Thus we chose to undertake the study using a spatial resolution of 1 km 2 .
We used pseudo-absences to use our species distribution models. We selected pseudo-absences following three criterion concerning the total number of pseudo-absences, the number of locations of pseudo-absences from a given cell of 1 km 2 , and the location of pseudo-absences in forest areas (Sect. S1 in the Supplement).

Model selection for the lynx
We developed species distribution models based on the permanent occurrences data. The models were based on information-theoretic methods, which focus on the search for a parsimonious model as the primary philosophy of statistical inference (Burnham and Anderson, 2002;Johnson and Omland, 2004). We identified a set of a priori hypotheses on the estimated probability of the lynx occurrence that describe the natural conditions and resources required for refuge, food and breeding of the lynx based on knowledge on ecology and biology that exists for the Balkan lynx and Eurasian lynx (Table 1). We split the variables for species distribution modelling into two categories, namely "natural" and "human", e.g. a "two dimensional" habitat model that was developed for the brown bear in Spain (see Naves et al., 2003). The "natural" hypothesis has environmental variables that are related with food abundance of the lynx and forest connectivity. The "human disturbance" hypothesis assumed that human activities (e.g. land-use activities and roads) affect (lynx) mortality (Naves et al., 2003). Natural factors consisted of layers such as forest cover, elevation and the terrain ruggedness index (derived from slope of terrain), coniferous forests, broadleaved forests mixed forests, pastures, bare rocks and transitional shrub land-woodland composing model 1. The human variables comprised layers such as agricultural land, urban land, burnt land, Euclidean distance to road (i.e. proximity to roads), Euclidean distance to human settlements, and village density composing model 2 (see Sect. S2 for details). The layers depicting land cover, forest cover, roads and villages were transformed into a set of neighbourhood variables (calculating the mean value of the original variable within a specified neighbourhood radius around the target cell) from 1 to 15 grid cells (corresponding to a radii of 1, 2, 3, 4, 5, 10, 15 km) and resulted in a total of 84 neighbourhood variables for potential use in the lynx models. The neighbourhood variables helped obtain information about the scale at which the lynx perceived the landscape and at which (natural) resources needed to be available. For details on neighbourhood variables see Sect. S3. The range of neighbourhood values (i.e. 1, 2, 3, 4, 5, 10, 15 km) corresponded roughly to the known variation in home-range size of the lynx (Bojovič, 1978;Balkan Lynx Strategy Group, 2008). We removed all variables that were highly correlated (Pearson correlation test; r > 0.7) and variables that did not show statistically significant differences between occurrence and pseudo-absence locations of the lynx (Kruskal-Wallis test; p > 0.05) (Sect. S4). We calculated the spatial autocorrelation of the dependent variable and for further details on spatial autocorrelation see Sect. S5.

Model fitting
We used generalised linear models (GLM) with logit-link to relate occurrence and pseudo-absence data to sets of explanatory variables to predict the probability of occurrence of the lynx in the study area. GLMs are an extension of classic linear regression models (McCullagh and Nelder, 1989). We used binomial error structure (logistic regression) to account for the lynx data structure. A logistic regression model predicts the probability of the occurrence of the lynx (e.g. for food and or refuge) at a given location within the study area. An estimated high probability of occurrence of the lynx would indicate suitable habitat of the lynx (suitable habitat is above the threshold value of the probability of occurrence of the lynx that is estimated by the species distribution modelling; see Model selection for the lynx and Model evaluation). All GLM (and resource selection probability function (RSPF) models that calculate species habitat selection are shown in Model evaluation) were fitted within the program R 2.9.0 (R Development Core Team, 2009). In total, we fitted 11 GLM candidate models for the lynx comprising the two categories of variables (candidate models are in Sect. S6 and the results of best models are in Sects. S7, S8). We selected the most parsimonious model from 11 candidate models using the corrected Akaike information criterion (AICc) for small samples (Johnson and Omland, 2004). The uncertainty was assessed using AIC weights (e.g. Fer-nández et al., 2006). AIC weights represented the relative likelihood of a model. AIC weights was taken as the approximate probability that each model is the best model out of the set of all proposed models (Anderson et al., 2000). We also calculated the difference of corrected Akaike information criterion AICc between competing candidate models (in Sect. S7). We mapped out predictions of probability of lynx occurrences for the most parsimonious model to identify the high probability values of lynx occurrences (most suitable habitat) and the low probability values of lynx occurrences (marginal and non-suitable habitat). We calculated the average value of estimated probability values of the lynx occurrences from 0.50 (e.g. Naves et al., 2003) to 1.00 to select a threshold (e.g. Liu et al., 2005). An estimated probability value above the threshold value, e.g. of 0.5 would be a suitable habitat. Marginal habitat and non-suitable habitats had values below or equal to this threshold value, e.g. of 0.5 for the lynx. The layers of model 1 and of model 2 were used as the biodiversity layer and cost layer, respectively, in Zonation to identify conservation areas for the lynx.

Model evaluation
We evaluated our models using three methods. First, all models were evaluated by calculating the area under the receiver operating characteristic curve (AUC). The value of AUC is from 0 to 1. A model with a value of AUC ≥ 70 % performs well and a model with a high value of AUC ≥ 90 % indicates outstanding discrimination between estimated high probability values of lynx occurrences (occurrences) and estimated low probability values of lynx occurrences (pseudo-absences in this study) (Hosmer and Lemeshow, 2000). We also calculated the predictive accuracy of the most parsimonious (from all candidate models) model using the deviance explained in percentage (D 2 ). We used cross-validation for logistic regression with a binary dependent variable because we had only one data set to check if there was over-fitting in the fitted models (e.g. Fernández et al., 2003;Kanagaraj et al., 2011). Our data set was divided into 10 parts (folds) where nine folds were used for fitting the model and the tenth fold was used for testing the model. This was repeated 10 times so that each fold was used for model testing. We considered the value of predicted probability (i.e. "cross-validation estimate of accuracy") which was calculated as an average value of probability after testing the model 10 times. A predicted probability value is from 0 to 1 (for details see DAAG package in R Development Core Team, 2009). Second, we calculated the resource selection probability function (RSPF) Lele et al., 2013;Thurfjell et al., 2014) for 11 (GLM) candidate models. RSPF shows the probability that an available natural resource unit, e.g. chamois prey, is selected (Thurfjell et al., 2014) by the lynx when encountered. We assessed the model fit (logistic regression with 39 permanent occurrences and 39 pseudo-absences) for 11 RSPF candidate models using Hosmer-Lemeshow test for goodness of fit by calculating the log-likelihood, the Bayesian information criteria (BIC), and the chi-squared test statistic (χ 2 ) (Hosmer and Lemeshow, 2000). We calculated the p value for 11 RSPF candidate models to assess the statistical significance of 11 RSPF candidate models (see, e.g. . We ran the best RSPF model (out of 11 RSPF candidate models) with 10 permanent occurrences and 10 randomly selected pseudo-absences i.e. we removed (at random) 25 % of 39 permanent occurrences of the lynx. We calculated the Spearman correlation coefficient between the best RSPF model with 10 permanent occurrences and 10 randomly selected pseudo-absences and the best RSPF model with 39 permanent occurrences and 39 randomly selected pseudo-absences (i.e. no permanent occurrences removed). We removed (at random) 51 % of 39 permanent occurrences and of 39 pseudo-absences of the lynx. We ran the best RSPF model with 20 permanent occurrences and 20 randomly selected pseudo-absences of the lynx. The Spearman correlation coefficient was then calculated between this RSPF model with 20 permanent occurrences and 20 randomly selected pseudo-absences of the lynx and the best RSPF model with 39 permanent occurrences and 39 randomly selected pseudo-absences of the lynx. A high value of Spearman correlation coefficient shows a strong prediction by the best RSPF candidate model (e.g. . We then applied the best RSPF model in the Step-Selection Function (SSF) using Geospatial Modelling Environment (GME) (www.spatialecology.com/gme/). The SSF is normally used with telemetry data of species involving randomly walking of animal (see, e.g. Proulx et al., 2013). The SSF links consecutive 39 lynx locations that are called (used) steps (Thurfjell et al., 2014). We made three assumptions for the SSF. Firstly, we assumed that our lynx data were as accurate as the telemetry data of species for our data modelling (e.g. Johnson et al., 2008). Secondly, we assumed no seasonal variation for collected lynx data in year 2006. We expected five "available" random steps for each of the lynx's steps used (e.g. lynx canadensis in Thurfjell et al., 2014). For each used and random steps, we calculated the step length and turning angles using non-linear functions (for details, see Beyer, 2015). Thirdly, we assumed an exponential function to calculate the lynx step length with a different single rate parameter and we assumed wrapped Cauchy function (see Beyer, 2015) to calculate lynx turning angles. We simulated our SSF model 30 times. We calculated the number of lynx steps predicted by each SSF model simulation. We then overlapped the predicted lynx steps by each SSF model simulation with predicted probability of occurrence of the lynx (above the threshold of estimated probability of occurrence of the lynx; see Model fitting) by our best GLM model. We selected the SSF model simulation with the highest number of lynx steps overlapped with predicted probability of occurrence of the lynx (above the threshold of the estimated probability of occurrence of the lynx) to calculate the mean of step length of the lynx. We compared the step length mean (in km) of the selected SSF model simulation with the range (radius in km) of neighbourhood variables of our best GLM models. We also overlapped the predicted lynx steps by the selected SSF model simulation with the prioritised conservation areas of 10 and 20 % by Zonation. Third, we randomly removed 1, 2, 4, 8, and 16 permanent occurrences and replaced the same number of removed permanent occurrences (1, 2, 4, 8, and 16) with randomly selected locations (no presence) from 136 cells of 10 km × 10 km of the lynx survey of  and we ran the best GLM model. We did the same for combined permanent and temporal lynx occurrences (see Sect. S9 in the Supplement).
Finally, we collected information on evidence for lynx presence (photos) in the study area. We used Google Earth, ArcGIS and information about the locations evidence of lynx presence gathered in the study area to overlap with the predicted probability of occurrence of the lynx (see Evaluation of GLM models). This information was collected in National Park Rugova in , western Kosovo, in Munella Mountain in 2014, in Thethi National Park, and in Thira in 2015, northern Albania, in Pelister National Park in 2013, and in Mavrovo National Park in 2015, Macedonia (Balkan Lynx Recovery Programme, 2014. In Albania, all lynx photos in 2014 were taken by camera traps that were placed from 1 to 2 km from each other covering an area of 400 km 2 (Balkan Lynx Recovery Programme, 2014).

Prioritisation of landscape for lynx conservation
The conservation planning software package "Zonation" (Moilanen et al., 2012) was used to identify well-connected habitat in the landscape suitable for lynx conservation (see Sect. S10 for the Zonation configuration settings used). Although Zonation is designed for prioritising areas with multiple species, it is also appropriate to use with a single species. Here, we utilise features of Zonation which allow us explore locations for extending existing conservation areas while accounting for connectivity (of the lynx habitat). The main out-put of Zonation is a landscape prioritisation in the form of a raster map where each raster cell is ranked in terms of its importance in providing lynx habitat. Thus, Zonation can be used to identify locations in the landscape that should be protected to ensure the conservation of the lynx. Using Zonation we examined the extent to which existing protected areas in the four countries overlapped the most important areas for lynx conservation as predicted by Zonation. We then undertook a separate Zonation analysis to examine how these protected areas should be expanded to maximise conservation benefits to the lynx.
We used the species distribution model derived from natural factors as the primary input for Zonation. The core-area variant of the Zonation algorithm was used (Moilanen et al., 2005), which minimizes biological loss by iteratively picking a cell i, for removal that has the smallest biological loss, δ until all cells have been removed. For a single species, the biological loss is determined by The function Q i (S) identifies the proportion of the remaining distribution of lynx habitat located in cell i, out of in the remaining cells to be prioritised, S. When a part of the distribution of the species is removed by Zonation, the proportion located in each remaining cell is raised, which results in Zonation attempting to maintain core areas of lynx habitat until the process of cell removal was finished. C i is the cost associated cell i (Moilanen, 2007). We used the species distribution model derived from human factors to generate the cost layer. This was chosen as all the factors in this model were thought to be detrimental to the lynx, and we inverted the layer so that areas with low values of the estimated probability of the lynx occurrence (marginal and non-suitable habitat; see Model selection for the lynx and Model evaluation) had a high cost, and areas with high values of estimated probability of the lynx occurrence (suitable habitat) had a low cost of lynx conservation. (This means that suitable habitats are less expensive to maintain, manage, and conserve the lynx habitat located in natural and semi-natural forests compared to nonsuitable habitats located in agricultural areas for the lynx. We assume the land conversion from non-preferred agricultural land to preferred forested land would be associated with a high cost for the conservation of the lynx). The boundary quality penalty (BQP) function of Zonation was used in our analysis to incorporate spatial aggregation into the areas prioritised for the lynx for conservation. When using the BQP, Eq. (1) is replaced by a more complex equation to determine the cell removal rule (for details see Moilanen, 2007). This allowed the habitat value of each cell to be adjusted based on proportion of the cells containing habitat within a specified neighbourhood around a given focal cell. The BQP is characterised by an effect radius and a response curve. The effect radius specifies the size of the neighbourhood used in calculating the BQP. For the lynx we used an effect radius determined by the radius of the neighbourhood variable used in the best fitting model 1; this value was 4 km 2 ( corresponding to 50.2 km 2 ). The response curve specifies how the local quality of a focal cell changes when different proportions of areas with high probability of lynx occurrences are lost from within the effect radius, and it is related to the sensitivity of the lynx to fragmentation of its habitat (Moilanen and Wintle, 2007;Gordon et al., 2009). The BQP curve used for this analysis is given in the Supplement (Sect. S10). The Zonation analysis was done using maps with the same 1 km 2 resolution as used in the species distribution modelling.
Finally, Zonation was run in both "constrained" and "unconstrained" mode. In the unconstrained mode, the algorithm can pick cells purely based on Eq. (1). In the constrained model we included the current protected areas using Zonation's mask function (Moilanen et al., 2012). This forced the raster cells in the protected areas to have the highest rankings to determine areas to extend the current set of protected areas.

Results
The spatial autocorrelation at the model grain of 1 km for the occurrences of the lynx was statistically insignificant (p value < 0.001). The model 1 showed the highest accuracy (AUC = 82 %) and performed better than model 2 (AUC = 65 %) in terms of model selection (the lowest value of AICc) and accuracy ( Table 2). The higher value of AUC for the model 1 suggests that model 1 shows a relatively clear pattern of estimated high probability of the lynx occurrences with respect to natural factors (elevation, terrain index, chamois prey occurrences). The lower value of AUC for the model 2 shows that the estimated probability of lynx occurrences is more conditioned by natural factors (e.g. the abundance of chamois prey) than by human disturbance i.e. agricultural land. This model 1 has an AIC weight of 50 % showing that this model is likely the best performing model from the 11 GLM models (see Model selection of the lynx). The models performed best at the neighbourhood scale of 50.2 km 2 (a radius of 4 km). The model 1 neighbourhood scale (a radius of 4 km) implied that the lynx selected habitat (primarily for food and refuge), if the natural resources and prey were available at the range of approximately 50 km 2 . This scale matches to the estimations of home-range of the lynx (e.g. Bojovič 1978). This model 1 showed that the estimated probability of occurrence for the lynx increased with elevation (ME), steep terrain (TRIR), and areas with abundant prey (chamois, CHP) and decreased with wellconnected transitional land (TRAN) that presents the change of forest management from woodland (taller trees, higher tree volume) to shrub-land (shorter trees, lower tree volume) ( Table 2). Table 2. The summary of the GLM models of estimated probability of lynx occurrences used for the Zonation analysis and the model selection estimators. Model 1 composed of natural factors was used as the habitat layer in the Zonation analysis the output of model 2 composed of anthropological factors was inverted and used as the cost layer in Zonation. Neighbourhood area is the scale; AICc = corrected Akaike's information criterion; Akaike w i = Akaike weight; D 2 = deviance explained; CV = cross validation; AUC = area under curve i. For each model the coefficients and their sign is shown along with their standard errors and statistically significance.

Model
Neighbourhood  The selected model 2 indicated that the estimated probability of the occurrence for the lynx decreased in human modified landscapes (i.e. in agricultural land). Agriculture land had a negative impact and the greatest effect (Table 2) on the lynx habitat. The selected neighbourhood scale of this model 2 was 314 km 2 (a radius of 10 km). This scale indicated the scale at which the lynx avoids agricultural areas.
The models with combined permanent and temporal occurrences indicated that the lynx selected areas wellconnected with forests at the range of approximately 50 km 2 and avoid agricultural and urban areas at the range of 707 km 2 . The results of the best models with combined permanent and temporal occurrences are in Sects. S7-S9 and S11-S12 in the Supplement.

Evaluation of GLM models
The best RSPF model (selected model) was model 1 at the range of approximately 50 km 2 . This model fully agreed with the best GLM models in terms of sign of coefficients and of variables. RSPF model 1 is statistically significant (p value of 0.026) with a negative value of log-likelihood of 93.3, with a value of χ 2 of 17.3 and the lowest value of AIC of 196.6 from 11 RSPF models. The neighbourhood scale of this RSPF model 1 has a radius of 4 km, implying that the natural resources and prey are available and encountered by the lynx at the range of approximately 50 km 2 , the lynx selects these natural resources, e.g. for food (Table 3). The Spearman correlation coefficient between the best RSPF model (no permanent presence removed) and the best RSPF model, removing 25 % of the permanent occurrences, was 0.93. The Spearman correlation coefficient between the best RSPF model (no permanent presence removed) and the best RSPF model, removing 51 % of the permanent occurrences, was also 0.93, indicating strong prediction by the best RSPF model (see Model evaluation). We simulated SSF with the selected RSPF model (model 1). We selected SSF model simulation with 222 predicted lynx steps as the best SSF model simulation. We overlapped 222 predicted lynx steps by the best SSF model simulation with predicted probability of occurrence of lynx with H >0.65 (threshold) of best GLM model (model 1 in Table 2). The selected SSF simulation had the highest number of steps of 36 % and a mean step length of 3.3 km for 222 predicted lynx steps by 30 SSF model simulations. The mean step length of 3.3 km corresponds closely to the radius of 4 km of our neighbourhood variables (at range of approximately 50 km 2 ) for our best GLM models (Table 2). We mapped out the predicted probability of occurrence of the lynx (H > 0.65) by the best RSPF model and predicted probability of occurrence of lynx with H > 0.65 by our best GLM model. We found that maps of the best RSPF model and best GLM model were very alike (Fig. S9 in the Supplement). Predicted probability of occurrence of lynx with H > 0.65 by best RSPF model and predicted probability of occurrence of lynx with H > 0.65 by our best GLM model resulted in 11.9 and 11.8 % of the entire study area, respectively. This indicates GLM model adequately predicted the probability of occurrences of the lynx.
Finally, we found that predicted lynx steps by SSF model increased from 36 % in areas with high probability of lynx occurrence (H > 0.65) to 50 % in areas with low probability of lynx occurrence (H > 0.50 to H < 0.65), estimated by our best GLM model (Fig. 2c). We overlapped 222 predicted steps by the selected SSF simulation with prioritised conservation areas of 10 and 20 % by Zonation. We found approximately 20 and 56 % of all steps inside of prioritised conservation areas of 10 and 20 % by Zonation, respectively (Fig. 3). We also simulated SSF model with RSPF model 2 including agricultural land at the range of 314 km 2 . Then, we overlapped 195 lynx steps predicted by the SSF model (with RSPF model 2) simulation with predicted probability of occurrences of the lynx with H > 0.50 (threshold of estimated probability of occurrences of the lynx) for the GLM model 2 (see Table 2). We found that 89 % of all steps predicted by this SSF model (with RSPF model of anthropological factors) simulation were associated with no agricultural land.
The change of AICc was below 2 ( AICc <2) between models of the best GLM model (model 1 in Table 2) with 1, 2 randomly removed permanent occurrences and the best GLM model (no permanent occurrences removed). Models with randomly removed 4, 8, and 16 permanent occurrences had a change of AICc above 2 ( AICc > 2) compared to our best GLM model (see Model evaluation). Results are shown in Sect. S9.
All evidence of the lynx overlaps with areas with high probability of lynx occurrences of H > 0.65 and with areas with low probability of lynx occurrences (marginal habitat) of H > 0.50 and H < 0.65 in Fig. 2 and in Sect. S11. Four pieces of evidence numbered (1), (2), (5) and (6) were inside of protected areas, and numbers (3) and (4) were in the surrounding protected areas. Evidence piece (3) was 12 km, and evidence pieces (4) was a 5 km (straight line) far from the closest protected area boundary (see Fig. 2b).

Lynx habitat maps
We mapped out the estimated probability of occurrence of the lynx in the model 1 and model 2 (Table 2; Fig. 2). The cut-off value in the estimated probability of occurrence of the lynx (predicted habitat suitability, H ) in the model 1 was H > 0.65 and in model 2 H > 0.50. This cut-off enabled us to distinguish areas with high probability of lynx occurrence (e.g. for model 1 H > 0.65) from areas with low probability of lynx occurrence (e.g. for model 1 H < 0.65) (see, e.g. Naves et al., 2003). Areas with high probability of lynx occurrence (suitable habitat) tended to occur in forested areas and areas with high elevation, while areas with low probability of lynx occurrence (marginal and non-suitable habitat) was concentrated in low land and non-forested land (i.e. agricultural land; the expansion of agricultural land may decrease the areas of areas with estimated high probability of lynx occurrence of the lynx in the study area). Patches with estimated high probability of lynx occurrence were fragmented and divided by areas of estimated low probability of lynx occurrence.
The areas with estimated high probability of lynx occurrence was spatially distributed over all four countries in the study area. Two large areas with estimated high probability of lynx occurrence, one spanning Montenegro, Kosovo, and Albania, the other spanning Albania, Kosovo, and the FYR of Macedonia is a key finding of this analysis (indicated by the zoomed-in maps in (a) and (b) in Fig. 2). The areas with estimated high probability of lynx occurrence outside protected areas were identified within Albania, which may ensure the connection of two multi-border areas of Albania, Montenegro, and Kosovo in the north and Albania, Macedonia, and Kosovo in centre of the study area, respectively. Areas with estimated high probability of lynx occurrence were also identified in the cross-border between Macedonia and Greece.
The areas with probability of the lynx occurrence (H > 0.65) predicted by the model 1 and model 2, respectively, overlapped (i.e. map (c) and map (d) in Fig.2). The model 1 and model 2 showed insignificant value of correlation (Spearman correlation; r = 0.013). Areas of overlap in both models indicated the presence of forested elevated land (model 1) as well as less agricultural land (model 2). These areas tended to occur predominantly in Macedonia and Albania and in western Kosovo.
The areas with estimated high probability of lynx occurrence for the lynx were identified in protected areas that had forested elevated areas and species of lynx prey (model 1) and have little or no agriculture land (model 2). This indicated that protected areas were well-located for lynx conservation. The areas with estimated probability values from H > 0.50 to H < 0.65 in protected areas can still be beneficial for the lynx (e.g. the photos of lynx were collected in Pelister National Park; (6) in Map (b) in Fig. 2 in 2013 (Balkan Lynx Recovery Programme, 2014).

Prioritisation of the landscape for lynx conservation
Areas prioritised by Zonation (unconstrained solution) were concentrated in the centre and north of the study area with the majority of the top 20 % solution occurring in Constraining the Zonation solution with the current protected areas forced all protected areas to be in the top 20 % solution (and most in the top 10 %) (Fig. 3c). The constrained (Fig. 3c) and unconstrained solutions (Fig. 3d) were similar. This indicated again that the existing protected areas occurred in high priority locations for lynx conservation. The major differences between constrained and unconstrained solutions were a reduction in some of high priority areas in Albania, Macedonia, and between Albania and Macedonia border areas.
In the constrained solution, high priority rankings have been assigned to new areas outside the current protected areas in western Macedonia, the north and east of Albania as well as in northern and western Montenegro, (Fig. 3c). These new areas within the top 20 % landscape prioritisation connected the two protected areas situated in the Albania-Montenegro-Kosovo cross-border area (Fig. 3a) and the Albania-Macedonia-Kosovo cross-border area (Fig. 3b). Prioritised areas in Macedonia and Albania tended to be fragmented.
The number of steps predicted by SSF in top 10 % Zonation solution was larger for Macedonia compared to Albania, with Kosovo having a small share and no share for Montenegro (Fig. 4, light grey components of columns). The top 20 % solution increased the difference between Macedonia and three other countries (Albania, Montenegro, and Kosovo), with Macedonia having the largest share of 68.2 % (Fig. 4).

Discussion
We used a GLM with an information-theoretic approach for testing a set of a priory hypothesis on the estimated probability of occurrences of lynx. This modelling approach helped us to identify the factors that determine the distribution of the lynx. Our models with permanent occurrences (the number of positive answers was above 50 % affirmative answers collected on lynx by interviewers for a cell of 10 km × 10 km) and combined permanent and temporal lynx occurrences (the number of positive answers was below 50 % affirmative answers collected on lynx by interviewers for a cell of 10 km × 10 km) showed that distribution of the lynx is determined by the abundance of environmental resources, high elevation, slope, the presence of forested and pastureland, and of lynx prey (chamois). Other studies on the Iberian lynx (Lynx pardinus) have also found that the lynx in Spain was conditionally distributed by natural resources, shrubs and rocky areas (Fernández et al., 2006) and pastureland (Fernández et al., 2003). The GLM models are used in habitat modelling of large-carnivore species because they are an effective modelling approach for SDMs performing well with adequate species data and handling non-linear responses between species occurrences and environmental variables (Kanagaraj et al., 2013). There are many approaches of species distribution modelling such as Maxent (Phillips et al., 2006) that uses presence-only data of species, niche models like ecological niche factor analysis (ENFA) (Hirzel et al., 2004) that estimates species home-range, resource selection functions (RSFs) and resource selection probability functions (RSPFs) (Lele et al., 2013) that calculates habitat selection of species. GLMs deal with autocorrelation of species data, multi-collinearity of variables and with model performance using very few species occurrence records, though GLMs do not show disadvantages (see, e.g. Kanagaraj et al., 2013). Another issue concerning species distribution modelling including GLMs is on the concept of "habitat" (a good discussion is provided by Gaillard et al., 2010). In this study, we applied RSPF model (model 1) to calculate the probability of occurrence of the lynx to compare explicitly with probability of occurrence of the lynx estimated by our best GLM model (model 1). We found that results obtained by best RSPF model and by best GLM model were very similar in the terms of model type, i.e. "natural" vs. "human" model, number and type of explanatory variables and sign of coefficients (see Table 2) and the maps ( Fig. 2c and Fig. S9). These findings indicated that the predictions of RSPF model and GLM model were consistent with known information of lynx occurrences in the cross-border area between Albania and Macedonia  and of lynx occurrences in Albania, Macedonia, and Kosovo (see evidence of lynx collected from 2013 to 2015 in Albania and Macedonia shown in Figs. 2a, b and S11a, b).
We checked our lynx data (see Lynx data) quality and detectability of our best GLM model by randomly removing a different number of permanent occurrences and a different number of combined permanent and temporal occurrences and replaced removed lynx occurrences with no lynx occurrence records (see Model evaluation). We found that our best model performed consistently worse than our best GLM model ( AICc > 1; CV < 0.71) after we randomly removed permanent lynx occurrences (Sect. S9) showing that lynx permanent occurrences could be better identified than lynx temporal occurrences. Overall, the lynx showed clear patterns of areas with high probability of its occurrences with respect to natural resources. Our findings demonstrate that using fine resolution environmental data (elevation, forest, land) helps detect areas with high probability of occurrence of species even if we use coarse resolution species data. Our results of Step Selection Function (SSF) model showed that lynx prefers to inhabit mostly in the cross-border area of Albania, Macedonia, and Kosovo. We suggest SSF or any resource selection models that tackle the autocorrelation problems of telemetry data (see, e.g. Johnson et al., 2008) to be used with new radio-telemetry data of lynx that are collected (personal communication with Gjorge Ivanov) to understand the animal behaviour and movement. This could be a new research area for organisation teams (see Sect. S0) working on the lynx data collection and lynx conservation.
This study clearly demonstrated the conservation areas for the lynx. We assessed potential key conservation areas for lynx conservation in south-eastern Europe and provided a large-scale study of the estimation of probability of occurrence of the lynx in four post-socialist countries of Albania, Macedonia, Montenegro, and Kosovo. This is an important step forward for designing conservation plans for large carnivores in countries that experienced massive political and socioeconomic changes (Radeloff et al., 2013). Our study identified the two cross-border areas between Albania and Macedonia and between Albania, Montenegro, and Kosovo in Fig. 3a and b for conservation of this enigmatic species. These results suggest that a close and effective cooperation of the governments of Albania, Macedonia, Montenegro, and Kosovo is required for lynx conservation in these cross-border areas (shown in box (a) and (b) in Fig. 3). The cross-border area between Albania and Macedonia is the key-conservation area identified by local and international organisation in terms of forest conservation, of local population awareness and engagement in nature conservation and of nature conservation at national and international level (see Sect. S0.1).
Our model results using the lynx occurrence records from questionnaires are also supported by direct observation records of the lynx subspecies (box (a) and (b) in Fig. 2). Our results may also be used for identifying new areas for intensive research. This may also help meeting targets of the Convention on Biodiversity Conservation of 17 % of terrestrial ecosystems and of no unknown species to go extinction (e.g. the lynx) by 2020 by governments. These findings support well the cross-border cooperation that has already started between Albania and Macedonia in 1990s (see Sect. S0.1).
This paper showed a full workflow starting with lynx occurrence data, applying species distribution modelling and combining it with spatial prioritisation of locations for lynx conservation. Our study is added to the existing literature of conservation science that has used generalised linear models and Zonation (see, e.g. Kujala et al., 2013). We separated natural and human factors to estimate the probability of the lynx occurrences and then to map "natural" and "human" model (the lynx habitat maps were based on "natural" and "human" factors, respectively). These "natural" and "human" maps served as the biological feature layer and the cost layer ("human" model was inverted), respectively, in Zonation software. We assessed, thus, the conservation areas of the lynx. This helped to show where to extend the existing protected areas and to potentially create new protected areas (and buffer zones) for lynx conservation.
The needs of lynx for forests and pastures (pastures, prey, well-connected forests were potentially selected habitats of the lynx for refuge and food) could be jeopardised by humans. We observed disturbed forests near stable forests and in protected areas and their surroundings (Sect. S13). Similar results were found in the case of Ukraine and Romania (postsocialist countries). In Romania, the high logging rates were likely trigged by rapid changes in institutions and ownership that resulted in forest decrease inside protected areas, which in turn, caused an increasing fragmentation of forest cover in the protected areas (Knorn et al., 2012).

Recommendations for lynx conservation
This study showed lynx conservation areas and predicted lynx steps by SSF were spread between Albania and Macedonia (Figs. 3 and 4). This indicates that these countries have similar responsibilities for the conservation of the lynx populations and habitat. Albania and Macedonia currently have the largest lynx populations and both countries share the largest part of predicted steps of lynx by SSP in the top 10 % of the landscape prioritised for lynx conservation (Fig. 4). Because of the existing lynx populations in Albania and Macedonia, both countries will be required to make all efforts to maintain these populations. Albania and Macedonia can collaborate to ensure a single long-term conservation plan of their cross-border areas, for example, to collaborate on preparing a species conservation plan. In addition, the results of our species distribution modelling and landscape prioritisation indicate that both countries may need to address the following factors: (i) the conserving of existing areas with high probability of lynx occurrence through forest conser-vation measures such as incentive programs to reduce the forested areas for forest harvesting and firewood (e.g. firewood is still important source for heating in rural areas in Albania; Laze, 2014); such an incentive program may compensate local people for the switching from wood harvesting to forest conservation (i.e. from resource to biodiversity forest ecosystem service); (ii) well-managed hunting activities such as limiting or halting chamois, brown hare and roe deer hunting at fixed times of the year; (iii) sustainable planning for the use of agriculture and pasture land; (iv) the collecting of constant seasonal and yearly lynx data using the advance monitoring system (radio-telemetry); (v) expansion of protected areas and the addition of new protected areas in locations such as the cross-border areas but also inside of Albania, Macedonia, Kosovo, and Montenegro (Fig. 3). This last point could help the lynx disperse from the north to the centre of the study area. This is because any lynx individual matters for the existing small population of the lynx. (Data collection on Carpathian lynx occurrence in the north of the study area (primarily Montenegro) can be useful to identify any potential overlapped areas of the distribution between Carpathian lynx and the lynx; yet, this is not our study objective).
Our results suggest that conservation efforts in the countries of our study could be oriented towards increasing the size of existed protected areas and their surrounding-buffer areas. This, however, will require a closer inter-country collaboration because the most promising areas for lynx conservation extend over the country borders. The expansion of protected areas must be associated with the conservation of forests, the reduction of human disturbance in the prioritised areas as well as with the connectivity of habitat patches for the lynx in their cross-border areas. This is because the lynx prefers larger areas of undisturbed forests, moves a great deal within its home-range (Balkan Lynx Strategy Group, 2008) to satisfy its needs for refuge, food and breeding as well as it is sensitive to anthropogenic land changes. We highlight here the relevance of large-area for the production of the lynx quoting as follows: Although isolation is of great importance in the production of new species, on the whole I am inclined to believe that largeness of area is still more important, especially for the production of species which shall prove capable of enduring for a long period, and of spreading widely (Darwin, 2011).
The Supplement related to this article is available online at doi:10.5194/we-16-17-2016-supplement.