337 Identifying habitat use of Ursus arctos, Lynx lynx martinoi and Canis lupus lupus in Albanian forests using occupancy modelling Abstract Forests are the principal terrestrial ecosystem for protected and endangered large carnivores, providing the main habitat for species maintenance and survival. Changes in forest cover influence species distribution. The aim was (1) to test hypotheses on the natural environmental (abiotic) and biological (biotic) factors and human disturbances that determine the colonization and local extinction of three large carnivore species in relation to forest cover, (2) to infer the biotic interactions between these three large carnivore species occupying the same forested areas in Albania. Colonization is estimated to be positively affected by forest cover for brown bear, Balkan lynx and grey wolf. Brown bear and grey wolf tend to compete for the same forested areas. Local extinction increased with decreasing forest cover for brown bear and increased with decreasing mixed broadleaved forests for Balkan lynx. Anthropological variables (proximity to villages and to neighbourhood roads) increased local extinction for brown bear, grey wolf and Balkan lynx. Further studies are recommended for better understanding biotic interactions of large carnivore species in forest habitats in Albania and its neighbouring countries, which could contribute to conservation of large carnivore species on a large scale. Iz vleček Gozdovi so glavni kopenski ekosistemi za zavarovane in ogrožene velike zveri in predstavljajo osnovni habitat za ohranjanje in preživetje teh vrst. Pokrovnost z gozdom vpliva na razširjenost vrst. Namen raziskave je bil (1) testirati hipoteze o vplivu naravnih okoljskih (abiotskih), bioloških (biotskih) dejavnikov in motenj zaradi človeka, ki vplivajo na naseljevanje in lokalno izumiranje treh velikih zveri v odvisnosti od pokrovnosti gozda, (2) sklepati na biotske odnose med tremi vrstami zveri, ki naseljujejo ista gozdna območja v Albaniji. Večja pokritost z gozdom pomeni večje naseljevanje rjavega medveda, balkanskega risa in sivega volka. Medved in volk običajno iščeta za ista gozdna območja. Lokalno izumiranje medveda se je povečalo zaradi zmanjšanja gozdnih površin, lokalno izumiranje risa pa je se povečalo na račun zmanjšanja površine mešanih gozdov. Zaradi dejavnikov, odvisnih od človeka (bližina vasi in prometnic), se je povečalo lokalno izumiranje medveda, volka in risa. Za razumevanje biotskih odnosov med tremi vrstami velikih zveri v gozdnih habitatih v Albaniji in sosednjih državah so potrebne nadaljnje raziskave, kar bo omogočilo varovanje velikih zveri v večjem merilu. Key words: Balkan lynx, biotic interactions, brown bear, colonization, local extinction, grey wolf. Ključne besede: balkanski ris, biotski vplivi, rjavi medved, naselitev, lokalno izumrtje, sivi volk. Received: 21. 3. 2019 Revision received: 13. 2. 2020 Accepted: 20. 4. 2020 1 Polytechnic University of Tirana, Faculty of Civil Engineering, Department of Environmental Engineering, Rr. “M. Gjollesha”, No. 54, 1023 Tirana, Albania. E-mail: kuendalaze@gmail.com Kuenda Laze1  DOI: 10.2478/hacq-2020-0007 19/2 • 2020, 337–347 19/2 • 2020, 337–347 338 Kuenda Laze Identifying habitat use of Ursus arctos, Lynx lynx martinoi and Canis lupus lupus in Albanian forests using occupancy modelling Introduction Free populations of large carnivore species are valuable for conservation. Conservation of such populations requires adequate knowledge of the distribution of large carnivore species, of abiotic factors affecting their populations and of biotic factors influencing the coexistence or competi- tion of such species (e.g., Ordiz et al. 2015; Heim et al. 2017) on a large scale. Species distribution models (SDMs) provide informa- tion on how animals or plants selectively use natural re- sources. Natural resources selected by animals or plants are then incorporated as the main environmental factors into species models. For example, Maxent (Phillips et al. 2006) relies on the assumption that a focal species is de- tected everywhere the species is present (Yackulic et al. 2013) and no information on species absences is avail- able. However, a species cannot always be detected at a given site during a data collection survey. A species may not be present at these sites (Araújo & Guisan 2006) or the species may be overlooked (Kéry 2011) during field- work. Missing information about these sites may result in underestimated species distributions (see, e.g., Kéry & Schaub (2012)), representing a gap in SDMs. Dynamic occupancy modelling fills this gap in SDMs by consider- ing various time periods and sites of animal distribution. When dealing with free species populations, species de- tectability of less than one leads to false negatives, where- by animal species are present but not detected during the survey (Kéry 2011). Species distribution models assume that a species selects the most preferred resources, though there are limitations, e.g., small species populations, dis- tant resources or human presences that prevent species from selecting the preferred resources (see, e.g., Araújo & Guisan (2006)). Dynamic occupancy models (and oc- cupancy models in general) make it possible specifically to account for the effect of local processes of coloniza- tion and local extinction on the species range dynamics (MacKenzie et al. 2003). Dynamic occupancy models have been successfully used to identify environmental or anthropogenic drivers in habitat use by brown bear in the French Pyrenees (Piédallu et al. 2019), to explain wolf recolonization in France (Louvrier et al. 2017) and to evaluate Eurasian lynx populations in the European Alps (Molinari-Jobin et al. 2017). Three protected large carnivore species were investigated in this study – brown bear (Ursus arctos), the critically en- dangered Balkan lynx (IUCN 2015) (Lynx lynx martinoi), which is a subspecies of Eurasian lynx, and grey wolf (Canis lupus lupus) in Albania. Populations of brown bear, Balkan lynx and grey wolf are protected by the national law “On biodiversity protection” of 2006 and the National Red List of Flora and Fauna of Albania (nationalredlist.org), and the Bern Convention Appendix II (coe.int/en/web/bern-con- vention). Brown bear and grey wolf populations in Albania are important for connecting brown bear populations in Europe and grey wolf populations in the Balkan Peninsula (Kaczensky et al. 2012b). The Balkan lynx population is mainly found in the cross-border area between Albania and North Macedonia. Estimates of large carnivore species pop- ulation sizes slightly increased from approximately 180 to 200 brown bears, from 200 to 250 grey wolves from 2006 to 2010/2011 (Kaczensky et al. 2012b) and decreased from 25 to 5 Balkan lynxes between 2001 and 2011 in Albania (von Arx et al. 2004; Kaczensky et al. 2012b). Species occupancy was assessed by estimating the prob- ability of initial occupancy, extinction, recolonization and detection for our large carnivore species in the years 2001 and 2011. A set of drivers of colonization and local ex- tinction was tested and potential interspecific interactions were inferred for the three protected large carnivore spe- cies in forests in Albania between the years of 2001 and 2011. The most recent literature on brown bear, lynx and wolf in Europe and the best available online data of natu- ral environment (abiotic) and biological (biotic) factors were tested for brown bear, Balkan lynx and grey wolf in Albania. We used large carnivore species occurrence data from 2001, forest cover from 2000 (Hansen et al. 2013), land use data from 2001, large carnivore species data from 2011 (Kaczensky et al. 2012b) and forest data from 2010 (Broxton et al. 2014). Specifically, the aim was (1) to understand the natural environmental (abiotic), bio- logical (biotic) and human disturbance factors that de- termine the probability of colonization and local extinc- tion of large carnivore species in forests in Albania and (2) to make inferences about biotic interactions of large carnivore species occupying the same forested areas. This study may help provide insights into biological biotic fac- tors that influence the coexistence or competition of large carnivore species and potential abiotic factors influencing species population change (future expansion or reduc- tion) of brown bear, Balkan lynx and grey wolf in Al- bania, with possible implications in neighbour countries. Materials and Methods Study area and large carnivore species data The study area was the entire territory of Albania. Brown bear and grey wolf data were collected in Albania in 2001 by the EMERALD I and II projects (Council of Europe and Ministry of Environment Forests and Water Admin- 19/2 • 2020, 337–347 339 Kuenda Laze Identifying habitat use of Ursus arctos, Lynx lynx martinoi and Canis lupus lupus in Albanian forests using occupancy modelling istration 2006). Balkan lynx data were provided by Balkan Lynx Survey Europe 2001 – Balkan population conduct- ed by KORA (a non-governmental organisation based in Switzerland) (von Arx et al. 2004) and the Interna- tional Union for Nature Conservation (IUCN), IUCN/ SSC Cat Specialist Group (IUCN/SSC Cat Specialist Group 2012). Species data for brown bear, Balkan lynx and grey wolf from 2011 were prepared and published online by the IUCN/SSC Large Carnivore Initiative for Europe for the European Commission (Kaczensky et al. 2012a, b). The large carnivore species distribution maps were composed of ‘sporadic occurrences’ and ‘permanent occurrences’ with a resolution of 10 km × 10 km. Large carnivore species permanent occurrence records (hence- forth ‘occurrences’) for 2001 were considered for brown bear, Balkan lynx and grey wolf and permanent occur- rences records for 2010/2011 for brown bear and grey wolf and sporadic occurrences for 2011 for Balkan lynx. There were 24 permanent occurrences for Balkan lynx, 41 for brown bear and 93 for grey wolf for 2001 and 12, 66 and 125 occurrences records for Balkan lynx, brown bear and grey wolf (Figure 1), respectively, in 2011. Figure 1: Locations of permanent occurrences of brown bear, Balkan lynx and grey wolf in Albania. Large carnivore species data for 2001 were collected by KORA for the Eurasian Lynx Online Information System for Europe (von Arx et al. 2004). Large carnivore species data for 2011 were available from the European Commission (Kaczensky e tal. 2012b). Slika 1: Lokacije stalnega pojavljanja rjavega medveda, balkanskega risa in sivega volka v Albaniji. Podatki o velikih zvereh za leto 2001 so zbrani v KORA za Eurasian Lynx Online Information System for Europe (von Arx et al. 2004). Vir podatkov za leto 2011 je Evropska komisija (Kaczensky et al. 2012b). Occurrences records 2001 2011 Balkan lynx Occurrences records 2001 2011 Brown bear Tirana Tirana Occurrences records 2001 2011 Grey wolf Tirana 19/2 • 2020, 337–347 340 Kuenda Laze Identifying habitat use of Ursus arctos, Lynx lynx martinoi and Canis lupus lupus in Albanian forests using occupancy modelling Predictor variables Abbreviation Species Model Assumption Spatial resolution Land cover 2001 Beech pure and mixed with coniferous forests (percent, km 2 ) 1 Brown bear, Balkan lynx ψ 1 , γ, ε, p + 15 m Mixed broadleaved class (percent, km 2 ) 2 Brown bear, Balkan lynx ψ 1 , ε + 15 m Coniferous forests (percent, km 2 ) 3 Brown bear, grey wolf ψ 1 , γ, p + 15 m Cultivated land (percent, km 2 ) 4 Brown bear, Balkan lynx, grey wolf ε - 15 m Mediterranean macchia (percent, km 2 ) 5 Brown bear, Balkan lynx, grey wolf ψ 1 , ε + 15 m Bare rocks and soil (percent, km 2 ) 6 Balkan lynx ψ 1 , p + 15 m Forest cover Forest & woodland cover 2000, (percent, km 2 ) 7 Brown bear, Balkan lynx, grey wolf ψ 1 , ε, p + 30 m Forest & woodland cover 2010, (percent, km 2 ) 8 Brown bear, Balkan lynx, grey wolf γ, p + 245 m Natural environment Elevation (Meters) 9 Brown bear, Balkan lynx, grey wolf ψ 1 , γ, ε, p + 15 m Terrain ruggedness index, (Unitless) 10 Brown bear, Balkan lynx, grey wolf ψ 1 , γ, p + 15 m Biological Brown bear neighbour grid cells (10 km× 10 km) 11 Brown bear γ, p + 1 km 2 Balkan lynx neighbour grid cells (10 km× 10 km) 12 Balkan lynx γ, p + 1 km 2 Grey wolf neighbour grid cells (10 km× 10 km) 13 Grey wolf γ, p - 1 km 2 Anthropogenic Distance to nearest dwelling road, (Meters) 14 Brown bear, Balkan lynx, grey wolf ψ 1 , ε, p + 1 km 2 Distance to nearest village, (Meters) 15 Brown bear, Balkan lynx, grey wolf ψ 1 , ε, p + 1 km 2 Road density, (km) 16 Brown bear, Balkan lynx, grey wolf ψ 1 , ε, p - 1 km 2 Table 1: List of sixteen predictor variables with spatial resolution of raw data as shown in the ‘Spatial resolution’ column used for the occupancy models, after predictor variables had passed the Pearson correlation and Kruskal-Wallis tests (S2). The data resolu- tion of this study was 1 km 2 . Assumptions on relationships between species colonization probability and predictor variables are based on the literature and knowledge of biology for brown bear (Piédallu et al. 2019), lynx (Molinari-Jobin et al. 2017; Melovski et al. 2018) and wolf (Louvrier et al. 2017) in Europe (see 2.3 Model selection). ψ 1 = initial occupancy probability γ = colonization probability, ε = local extinction probability, and p = detection probability, Y = year, 7 = forest cover in 2000, 8 = forest cover in 2010. All predictor variables used for the occupancy models are shown in Table S1. Hypotheses are shown in Table 2. Tabela 1: Seznam šestnajstih spremenljivk s prostorsko resolucijo surovih podatkov prikazanih v stolpcu ‘Spatial resolution’, ki smo jih uporabili za modele zasedenosti prostora, po testiranju spremenljivk s Pearsonovim korelacijskim testom in Kruskal- -Wallisovim testom (S2). Ločljivost podatkov je 1 km 2 . Predvidevanje odnosov med verjetnostjo naselitve in spremenljivkami je ocenjena na osnovi literature in poznavanja biologije medveda (Piédallu et al. 2019), risa (Molinari-ZJobin et al. 2017; Melovski et al. 2018) in volka (Louvrier et al. 2017) v Evropi (glej poglavje 2.3 Model selection). (ψ 1 = začetna verjenost zasedenosti prostora γ=verjetnost naselitve, ε= verjetnost lokalnega izumiranja, in p=verjetnost detekcije, Y=leto, 7=površina gozda leta 2000, 8= površi- na gozda leta 2010. Vse spremenljivke v modelih zasedenosti prostora so prikazane v Tabeli S1. Hipoteze so prikazane v Tabeli 2. 19/2 • 2020, 337–347 341 Kuenda Laze Identifying habitat use of Ursus arctos, Lynx lynx martinoi and Canis lupus lupus in Albanian forests using occupancy modelling Twelve sporadic occurrences of Balkan lynx in Albania (a grid cell resolution of 10 km × 10 km that was less than 50 percent constantly occupied by lynx (von Arx et al. 2004)), were assumed to be permanent occurrences (a grid cell of resolution of 10 km × 10 km that was more than 50 percent constantly occupied by lynx (von Arx et al. 2004)) for Balkan lynx because they were connected to permanent cross-border Balkan lynx occurrences in North Macedonia (Kaczensky et al. 2012b). On the assumption that large carnivore species may be present anywhere within a 10 km × 10 km grid cell, and taking into account finer resolution environmental variables (see 2.2. Environmental variables), the large car- nivore species occurrence records (10 km × 10 km) were downscaled to a 1 km × 1 km grid cell within the original (10 km × 10 km) grid cell (employing the Spatial Join function of the Analysis T ool in ArcGIS). Large carnivore species presence 1 km × 1 km grid cells were randomly assigned for each 10 km × 10 km grid cell as a species occurrence (Figure 1). Environmental variables Data on topography, land cover in 2001 and in 2010, forest cover in 2000 and 2010, and proximity and density of villages and roads were compiled (Table 1). Cultivated land and Mediterranean macchia shrubland data were from the Albanian National Forest Inventory (Agrotec. SpA.Consortium 2004a) with a spatial resolution of 15 m for 2001. The cultivated land and Mediterranean macchia data set was derived from satellite images of Landsat 5 TM and Landsat 7 ETM+. The accuracy of the land use and land cover maps was 67 percent for thicket and shrub lands (Agrotec.SpA.Consortium 2004b). Technical details of land classifications are described in (Laze 2013). ‘Woody savannas’ was used as a proxy for Mediterra‑ nean macchia shrubland (for 2001) provided by Broxton et al. (2014) for 2011. The data were derived from 10 years (2001–2010) of MODIS-based Global Land Cover Climatology (Collection 5.1 MCD12Q1), with a spa- tial resolution of 245 m. For each 1 km 2 pixel, the cor- responding shrubland cover class with the highest overall confidence was selected. The percentage forest cover was used for 2000 with a resolution of approximately 900 m 2 provided by Hansen et al. (2013). Trees were defined as vegetation taller than 5 m height (Hansen et al. 2013) in the forest cover class. Elevation (15 m × 15 m) and terrain ruggedness index (TRI) at a resolution of 1 km² (Naves et al. 2003) were derived from AsterDEM, which was a product of METI and NASA (METI and NASA 2011). The vector lay- ers of roads and villages were provided by Environmen- tal Legislation and Planning Albania. Neighbourhood roads (dwelling roads in Table S1) and other road types (asphalted, well-kept, seasonal roads) were taken into ac- count. Roads and village data were transformed as density (roads) and distance variables. Distance neighbourhood roads and distance villages variables were respectively the minimum Euclidean distances of the 1 km 2 cell to the nearest neighbourhood roads type and to the nearest vil- lage. The density variables was the densities of all road types per km² grid cell. Model selection A dynamic occupancy model (MacKenzie et al. 2003) was used to account for colonization probabilities, γ, representing the probability of a large carnivore species occurrence that was not occupied in 2001 becoming oc- cupied in 2011 and for extinction, ε, representing the probability of a large carnivore species occurrence that was occupied in 2001 becoming unoccupied in 2011, the initial occupancy ψ, representing the probability of a large carnivore species occurrence being occupied in the first year i.e., 2001, and the detection probability p, which is the probability of detecting large carnivore species occurrences. Large carnivore species movements were assumed to be random in preferred forested areas and their land surroundings (see Table 2) in 2000/2001 and 2010, and large carnivore species data were collect- ed twice (in 2001 and 2011), excluding large carnivore species breeding seasons, large carnivore species data uncertainty and large carnivore species detection un- certainty issues. Species occupancy can be interpreted as habitat use (suitability), which is the proportion of (grid cell) area occupied by the (large carnivore) species (MacKenzie & Nichols 2004). In relation to dynamic occupancy models for brown bear, lynx and grey wolf in Europe, a set of natural conditions and resources (natural environment), and biological and human dis- turbance predictor variables were selected for brown bear, Balkan lynx and grey wolf in Albania (Table 1). Occupancy models were built-in for the three large car- nivore species, focusing on time-varying natural condi- tions and resources predictor variables representing forest cover for 2000 and 2010, and 2001 and 2011 for three large carnivore species, and site-varying variables. Site- varying predictor variables were elevation, terrain rugged- ness index, shrub cover (Piédallu et al. 2019) and brown bear neighbouring grid cells, for brown bear. For lynx, site-varying predictor variables were elevation, terrain ruggedness index (Melovski et al. 2018) and lynx neigh- bouring grid cells (Molinari-Jobin et al. 2017). Terrain 19/2 • 2020, 337–347 342 Kuenda Laze Identifying habitat use of Ursus arctos, Lynx lynx martinoi and Canis lupus lupus in Albanian forests using occupancy modelling ruggedness index, shrub cover (Louvrier et al. 2017) and wolf neighbouring grid cells were the site-varying predic- tor variables for wolf. Site-varying human disturbance predictor variables were distance to nearest neighbour- hood roads and cultivated land for the three large carni- vore species (Louvrier et al. 2017; Molinari-Jobin et al. 2017; Piédallu et al. 2019). Prior to occupancy modelling, multicollinearity among predictor variables was evaluated using the Pearson corre- lation test (r > 0.7; S2). Spatial autocorrelation (similarity between occurrences record data of large carnivore spe- cies) of the dependent variables was calculated using the software Geoda095i (geodacenter.asu.edu), (S3). There were five time-varying predictor variable (years and forest cover) models and six site-varying predictor variable mod- els. Models were considered with either no yearly effect on colonization and on local extinction or no forest cover effect on detection for the three large carnivore species (e.g., Piédallu et al. (2017). Terrain complexity, occupan- cy of neighbouring grid cells for brown bear, occupancy of neighbouring grid cells for Balkan lynx, occupancy of neighbouring grid cells for grey wolf in 2011 and forest cover for 2000 and 2010 were the predictor variables for occupancy modelling of the distributions of brown bear, Balkan lynx and grey wolf in Albanian forests. Predictor variables were standardized prior to occupancy model- ling. Five time-varying and six site-varying predictor vari- able models were separately run using Akaike Informa- tion Criterion (AIC). The six site-varying models were composed of forest cover and terrain ruggedness index and then of shrubland, (pure) beech mixed with conifer- ous forests, mixed broadleaved forest, coniferous forest, road density, distance to nearest village and distance to nearest neighbourhood roads. One and zero, respectively, were assigned for occurrences (detection) and for non-de- tection in 2001 and 2011 for the six site-varying models, for every large carnivore species. The six occupancy mod- els with site-varying predictor variables for the three large carnivore species are shown in Table S4. Each occupancy model (11 models per year/per species) was separately run for each year (2001, 2011) and for each large carnivore species using the “unmarked 0.12.2” package (Fiske & Chandler 2011) in R v386 3.4.2 (R Development Core Team 2017). Hypothesis Description References Natural environment and biological factors Abiotic The species require natural dense, high stable and naturally less disturbed forests (by humans) and high elevation to hide, breed and search for food. Natural dense, high, stable and naturally less disturbed forest area should be available to brown bear, lynx and wolf positively affecting species colonization. Brown bear rely on oak and beech forests for food and lynx on bare rocks for breeding. (Manly et al. 2002; Zimmermann et al. 2005; Balkan Lynx Strategy Group 2008; May et al. 2008; Louvrier et al. 2017; Molinari-Jobin et al. 2017; Piédallu et al. 2017; Melovski et al. 2018; Grilo et al. 2019) Biotic Wolf prefers forested areas not used by brown bears (fewer brown bear occurrences) and not used by other wolf pairs (fewer wolf occurrences). Lynx tends to establish its home range near other lynxes (more lynx occurrences) and away from brown bears (fewer brown bears occurrences) in forests. (Krofel et al. 2012; Ordiz et al. 2015; Ražen et al. 2016; Molinari-Jobin et al. 2017) Human disturbance The species do not prefer cultivated land, villages and roads (high human disturbance). Human disturbance should negatively affect species colonization. (Naves et al. 2003; Ordiz et al. 2013; Zimmermann et al. 2014; Louvrier et al. 2017; Molinari-Jobin et al. 2017; Piédallu et al. 2017) T able 2: Hypotheses based on the literature for brown bear, Balkan and Eurasian lynx and grey wolf and knowledge of their biology. Tabela 2: Hipoteze na osnovi literature o rjavem medvedu, balkanskem risu in sivem volku ter poznavanju njihove biologije. 19/2 • 2020, 337–347 343 Kuenda Laze Identifying habitat use of Ursus arctos, Lynx lynx martinoi and Canis lupus lupus in Albanian forests using occupancy modelling Results No year or no forest cover time-varying predictor vari- ables effects were identified on initial occupancy, coloni- zation, local extinction and detection probability for the three large carnivore species. Models with time-varying predictor variables performed better for brown bear and for Balkan lynx, but worse for grey wolf compared to the models with site-varying predictor variables. No. Species Rank Model with time/site-varying predictor variables AIC ΔAIC AIC wi , (%) 1 Brown bear 1 ψ 1 (.), γ(.), ε(.), p(.) 150.4 0.0 51.6 2 Balkan lynx 1 ψ 1 (.), γ(.), ε(.), p(.) 39.8 0.0 51.5 3 Grey wolf 1 ψ 1 (.), γ(.), ε(.), p(.) 305.5 0.0 51.5 4 Brown bear 1 ψ 1 (7, 9, 15), γ(8, 10, 11), ε(7, 8, 15), p(9, 10, 8) 162.2 0.0 72.0 5 Balkan lynx 1 ψ 1 (7, 10, 14), γ(8, 10), ε(7, 14, 2, 1), p(10, 8) 39.8 0.0 38.1 6 2 ψ 1 (7, 10), γ(8, 10, 12), ε(7, 14, 2, 1), p(10, 8) 39.8 0.0 38.0 7 Grey wolf 1 ψ 1 (7, 9, 10), γ(8, 9, 10), ε(7, 8, 9, 10, 15), p(.) 30.1 0.0 38.1 8 2 ψ 1 (7, 10, 3), γ(8, 10, 13), ε(8, 5, 16), p(10, 8) 30.1 0.0 38.0 Notes: 1 = beech pure and mixed with coniferous forests; 2 = mixed broadleaved forests; 3 = coniferous forests, 5 = Mediter‑ ranean macchia, 7 = forest cover in the year 2000; 8= forest cover in the year 2010; 9 = elevation; 10 = terrain ruggedness index; 11 = brown bear neighbouring grid cells (10 km × 10 km), 12 = Balkan lynx neighbouring grid cells (10 km × 10 km), 13 = grey wolf neighbouring grid cells (10 km × 10 km), 14 = distance to nearest dwelling road, 15 = distance to nearest village, 16 = road density. Table 3: Summary of best-fit models first ranked with time-varying predictor variables and with site-varying predictor variables first ranked by Akaike’s Information Criterion (with ΔAIC < 2) for brown bear, Balkan lynx and grey wolf. ψ 1 = nitial occupancy probability γ = colonization probability, ε = local extinction probability, and p = detection probability, 7 = forest cover 2000 and 8 = forest cover in 2010. Best-fit models with time-varying predictor variables are numbered from one to three and best-fit models with site-varying predictor variables from four to eight. Six occupancy models variables, sign and magnitude of estimated coef- ficients, standard errors and p‑ value for the best-fit (first ranked only) occupancy models with site-varying predictor variables are shown in Tables S4 and S5. The local extinction probabilities of brown bear, Balkan lynx and grey wolf, estimated by the best-fit models, are shown in Figure 2. Tabela 3: Povzetek najbolj ustreznih modelov, razporejenih po časovnih in prostorskih spremenljivkah, razporejenih po Akaike’s Information Criterion (z ΔAIC < 2) za medveda, risa in volka. ψ 1 = začetna verjenost zasedenosti prostora γ=verjetnost naselitve, ε= verjetnost lokalnega izumiranja, in p=verjetnost detekcije, Y = leto, 7 = površina gozda leta 2000, 8 = površina gozda leta 2010. Najustreznejši modeli s časovnimi spremenljivkami so oštevilčeni od ena do tri in s prostorskimi spremenljivkami od štiri do osem. Šest modelov zasedenosti prostora, spremenljivke, značilnost in velikost ocenjenih koeficientov, standardna napaka in p-vrednost za najbolj usterezen model (samo prvo razporejeni) s prostorskimi spremenljivkami, so prikazani v Tablelah S4 in S5. Verjetnost lokalnega izumrtja medveda, risa in volka je ocenjena z najbolj ustreznim modelom, prikazanim na Sliki 2. Spatial autocorrelation was statistically insignificant (below < 0.007 in absolute terms) for the occurrences of the three large carnivore species in 2001 at 1 km². Forest cover had a positive effect on initial occupancy, colonization and detection probability for brown bear. Terrain ruggedness index positively affected colonization probability and detection probability, and high elevation positively affected the initial occupancy probability for brown bear. Brown bear neighbouring grid cells (which could be occupied by other brown bears) were positively associated with colonization probability. The local extinc- tion probability of brown bears decreased with decreasing human disturbance (i.e., greater distance to the nearest village) (Table S5). The initial occupancy and detection probability of Bal- kan lynx increased with increasing forest cover and with increasing terrain ruggedness index (terrain roughness). Forest cover and lynx neighbouring grid cells positively affected the colonization probability of Balkan lynx. Proximity to neighbourhood roads had a positive effect on Balkan lynx local extinction probability (Table S5). There was a positive effect of forest cover on coloni - zation probability and a positive effect of proximity to villages and road density on local extinction probability for grey wolves. The effect of forest cover and terrain rug - gedness index on initial occupancy probability was not estimated for grey wolf. No statistically significant predic- tor variables were found for grey wolf (Table S5). 19/2 • 2020, 337–347 344 Kuenda Laze Identifying habitat use of Ursus arctos, Lynx lynx martinoi and Canis lupus lupus in Albanian forests using occupancy modelling Discussion The availability of brown bear, Balkan lynx and grey wolf occurrence records collected at two time points (in 2001 and 2011), sample units (sites), and adequate information on locations of large carnivore species occurrences and on potential large carnivore species habitat unsuitability (ab- sences) in Albania allowed the use of occupancy models. Forest cover and forest cover connectivity determined the spatial distribution in relation to the spatial scale of home-ranges for brown bear, Balkan lynx and grey wolf in Albania using SDMs (Laze 2013; Laze & Gordon 2016). Changes in forest cover between 2001 and 2011 negative- ly affected the distributions of all three large carnivores in Albania (Laze 2019). This study provides new insights into the local processes of colonization and extinction of brown bear, Balkan lynx and grey wolf, as shown in Figure 2. Increasing forest cover positively affected the colonization probability of brown bear and the detection probability of Balkan lynx and grey wolf. Forest cover also positively affected colonization of (Canis lupus) wolf in France (Louvrier et al. 2017) and (Eurasian) lynx de- tection in Switzerland (Molinari-Jobin et al. 2017). This study showed that anthropological variables (distance to nearest villages and distance to nearest neighbourhood roads) increased the local extinction probability of brown bear and grey wolf, and Balkan lynx as shown by Table 2. Figure 2: Estimated probabilities of local extinction for brown bear, Balkan lynx and grey wolf. Estimated probability values of local extinction close to one were approximated to one. Slika 2: Ocenjena verjetnost lokalnega izumrtja rjavega medveda, balkanskega risa in sivega volka. Vrednosti ocenjene verjetnosti lokalnega izumrtja blizu ena, so bile zaokrožene na ena. Grey wolf Extinction probability 0 1 Balkan lynx Extinction probability 0 1 Brown bear Extinction probability 0 1 Tirana Tirana Tirana 19/2 • 2020, 337–347 345 Kuenda Laze Identifying habitat use of Ursus arctos, Lynx lynx martinoi and Canis lupus lupus in Albanian forests using occupancy modelling Decreasing mixed broadleaved forest cover was found to increase the local extinction probability of Balkan lynx. One explanation may be the changes in mixed broadleaved forest cover from 2001 to 2010. Changes in mixed broadleaved forest cover may also be associated with decreasing forest resource availability, e.g., loss of trees and loss of large carnivore prey inhabiting forests. Here, frequent, seasonal and standardized monitoring of forest data (old forests, young forests, tree height) and of large carnivore prey could increase data availability, pro- viding abundant information about forest cover, changes in forest cover and changes in the number of large car- nivore prey, helping to enhance our knowledge of sites visited by one or more large carnivore species (to be de- tected during a survey) in less naturally disturbed forests or intact forests. Adequate and abundant telemetry data of large carni- vore species during nursing time need to be collected in forests in order to understand better the biotic interac- tions of the three large carnivore species. Adequate time series killing, poaching (illegal killing) and density data of the large carnivore species should also be gathered in Albania. Killing of large carnivores can affect their species populations. Liberg et al. (2012) showed that wolf poach- ing accounted for half of total wolf mortality and over two-thirds of total wolf poaching remained undetected, showing that the wolf population in Sweden would have been almost four times as large in 2009 without poach- ing during the previous 10 years. Swenson et al. (1995) showed that brown bear populations were saved in Swe- den, but not in Norway, last century, because the most effective measures for conserving brown bear were those that reduced or eliminated the economic incentive for people to kill them. Linnell et al. (2017) showed that the Bern Convention provided scope for permitting the kill- ing of wolves (to address human-wolf conflicts), but ‘this must be clearly justified and proportional to the conser- vation status of wolves’ in order not to jeopardize wolf population recovery. Conclusions The distribution of brown bear, Balkan lynx and grey wolf was conditioned by forest cover. Loss of forest cov- er tended to increase the local extinction probability of these protected animal species. Brown bear and grey wolf tend to compete for the same forested areas in Albania. Regular data collection and monitoring of large carni- vore species, large carnivore prey and changes in forests is recommended for increasing information and knowl- edge of biotic interactions among species occupying the However, three anthropological variables differently in- fluenced the local extinction of each large carnivore spe- cies (villages for brown bear and grey wolf and neighbour- hood roads for Balkan lynx). Potential biotic interactions between brown bears and wolves Interspecific interaction is often not easy to incorporate into species distribution models (Godsoe & Harmon 2012). Dormann et al. (2018) demonstrated that species distribution models mostly assume biotic interactions to be constant in space and time, use data that fail to show real co-occurrences at ‘ecologically relevant spatial and temporal scales’, and cannot always use high-quality spe- cies distribution and environmental data and information from experiments for studying large carnivore species with similar or contrasting habitat requirements better. The re- sults of this study can show possible large carnivore spe- cies interactions by comparing the spatial distributions of the three large carnivore species, starting with interactions between brown bear and grey wolf in the same forests (see e.g., Ordiz et al. (2015)), which showed similar distribu- tion patterns. The changes in forest cover (between 2000 and 2010) may have changed the availability of forests, leading to potential increased competition between brown bears and grey wolves for food and refuge because they tend to select the same forested areas. For example, elevated forested areas tended to be occupied by both brown bears and grey wolves. Elevated forests (of 1 km 2 ) and complex topography areas may be selected by (Canis lupus) wolf pairs for refuge and breeding (Grilo et al. 2019) far from other grey wolf pairs (source populations) (Ražen et al. 2016). It can be inferred that brown bears could use the same forested areas for food by displacing grey wolf (pairs) from their kills. Ordiz et al. (2015) found that brown bear can affect (Canis lupus) (grey) wolf (and Eurasian lynx) by displacing wolf (and Eurasian lynx) from more than 50 percent of their kills, especially during grey wolf (and Eur- asian lynx) nursing time. Krofel et al. (2012) found that 32 percent of Eurasian lynx prey remains and 15 percent of large preys (biomass) killed by Eurasian lynx was taken by brown bear, pushing Eurasian lynx to increase the kill rate by 23 percent (in forest habitats). Data issues There were differences in data collection for large carni - vore species, forest and land cover available for Albania. Forest data were collected approximately in 2000 and 2010, land data in 2001, and large carnivore species were collected approximately in 2001 and 2010/2011. 19/2 • 2020, 337–347 346 Kuenda Laze Identifying habitat use of Ursus arctos, Lynx lynx martinoi and Canis lupus lupus in Albanian forests using occupancy modelling same forested areas and for upgrading priorities for forest and large carnivore species conservation. New large scale research could increase our knowledge about of interspe- cific competition among our focal large carnivore species in relation to forest habitats in Albania and its neighbour- ing countries. Kuenda Laze , https://orcid.org/0000-0002-2896-9548 References Agrotec.SpA.Consortium 2004a: Albanian National Forest Inventory (ANFI): Analysis of the spatio-temporal and semantic aspects of land cover/use dynamics. Agrotec SpA Consortium, Rome. Agrotec. SpA. Consortium 2004b: Albanian National Forest Inventory (ANFI): Special study on grazing impact on wooded lands, including fuelwood consumption assessment. Agrotec SpA Consortium, Rome. Araújo, M. B. & Guisan, A. 2006: Five (or so) challenges for species distribution modelling. Journal of Biogeography 33: 1677–1688. https://doi.org/10.1111/j.1365-2699.2006.01584.x Balkan Lynx Strategy Group 2008: Balkan Lynx Strategy Group. Strategy for the Conservation of the Balkan Lynx in Macedonia and Albania. Peshtani, MK, 3–4 June 2008. www.catsg.org/balkanlynx Broxton, P . D., Zeng, X., Sulla-Menashe, D. & T roch, P . A. 2014: A Global Land Cover Climatology Using MODIS Data. Journal of Applied Meteorology and Climatology 53: 1593–1605. https://doi. org/10.1175/jamc-d-13-0270.1 Council of Europe, Ministry of Environment Forests and Water Administration 2006: Second Emerald Project in Albania-Final report Dormann, C. F ., Bobrowski, M., Dehling, D. M., Harris, D. J., Hartig, F ., Lischke, H., Moretti, M. D., Pagel, J., Pinkert, S., Schleuning, M., Schmidt, S. I., Sheppard, C. S., Steinbauer, M. J., Zeuss, D. & Kraan, C. 2018: Biotic interactions in species distribution modelling: 10 questions to guide interpretation and avoid false conclusions. Global Ecology and Biogeography 27: 1004–1016. https://doi.org/10.1111/geb.12759 Fiske, I. & Chandler, R. 2011: unmarked: An R Package for Fitting Hierarchical Models of Wildlife Occurrence and Abundance. Journal of Statistical Software 43 (10). https://doi.org/10.18637/jss.v043.i10 Godsoe, W. & Harmon, L. J. 2012: How do species interactions affect species distribution models? Ecography (Cop) 35: 811–820. https:// doi.org/10.1111/j.1600-0587.2011.07103.x Grilo, C., Lucas, P . M., Fernández-Gil, A., Seara, M., Costa, G., Roque, S., Rio-Maior, H., Nakamura, M., Alvares, F ., Petrucci- Fonseca, F . & Revilla, E. 2019: Refuge as major habitat driver for wolf presence in human-modified landscapes. Animal Conservation 22: 59–71. https://doi.org/10.1111/acv.12435. Hansen, M. C., Potapov, P . V., Moore, R., Hancher, M., T urubanova, S. A., Tyukavina, A., Thau, D., Stehman, S. V., Goetz, S. J., Loveland, T. R., Kommareddy, A., Egorov, A., Chini, l., Justice, C. O. & Townshend J. R. G. 2013: High-resolution global maps of 21st- century forest cover change. Science 342: 850–853. https://doi. org/10.1126/science.1244693ž Heim, N., Fisher, J. T., Clevenger, A., Paczkowski, J. & Volpe, J. 2017: Cumulative effects of climate and landscape change drive spatial distribution of Rocky Mountain wolverine (Gulo gulo L.). Ecology and Evolution 7: 8903–8914. https://doi.org/doi:10.1002/ece3.3337 IUCN/SSC Cat Specialist Group 2012: Red list assessment. Statement of the IUCN/SSC Cat SG. http://www.catsg.org/balkanlynx/20_blx- compendium/home/index_en.htm Kaczensky, P ., Chapron, G., Von Arx, M., Huber, D., Andrén, H. & Linnell, J. 2012a: Status, management and distribution of large carnivores – bear, lynx, wolf & wolverine – in Europe. Part 1 Kaczensky, P ., Chapron, G., Von Arx, M., Huber, D., Andrén, H. & Linnell, J. D. C. 2012b: Status , management and distribution of large carnivores – bear, lynx, wolf & wolverine – in Europe Part 2 Kéry, M. 2011: Towards the modelling of true species distributions. Journal of Biogeography 38: 617–618. https://doi.org/10.1111/ j.1365-2699.2011.02487.x Kéry, M. & Schaub, M. 2012: Chapter 1 – Introduction. In: Kéry M, Schaub M (eds) Bayesian Population Analysis using WinBUGS. Academic Press, Boston Krofel, M., Kos, I. & Jerina, K. 2012: The noble cats and the big bad scavengers: effects of dominant scavengers on solitary predators. Behavioral Ecology and Sociobiology 66: 1297–1304. https://doi. org/10.1007/s00265-012-1384-6 Laze, K., 2013: Identifying and understanding the forest cover change patterns and processes in Albania and Kosovo. Halle, Univ., Naturwissenschaftlichen Fakultät III, Diss., 2013. Halle, Saale: Universitäts- und Landesbibliothek Sachsen-Anhalt. Germany Laze, K. 2019: Distribution of large carnivores (Ursus arctos, Lynx lynx, Canis lupus) in relation to spatial and temporal changes in forests: the case of Albania (Carnivora). Lynx new Ser 50: 61–86 Laze, K. & Gordon, A. 2016: Incorporating natural and human factors in habitat modelling and spatial prioritisation for the Lynx lynx marti- noi. Web Ecology 16: 17–31. https://doi.org/10.5194/we-16-17-2016 Liberg, O., Chapron, G., Wabakken, P ., Pedersen, H. C., Hobbs, N. T. & Sand, H. 2012: Shoot, shovel and shut up: cryptic poaching slows restoration of a large carnivore in Europe. Proceedings of the Royal Society B: Biological Sciences 279: 910–915. https://doi.org/10.1098/ rspb.2011.1275 Linnell, J. D. C., T rouwborst, A. & Fleurke, F . M. 2017: When is it acceptable to kill a strictly protected carnivore? Exploring the legal constraints on wildlife management within Europe’s Bern Convention. Nature Conservation 21: 129–157 Louvrier, J., Duchamp, C., Lauret, V., Marboutin, E., Cubaynes, S., Choquet, R., Miquel, C. & Gimenez, O. 2017: Mapping and explaining wolf recolonization in France using dynamic occupancy models and opportunistic data. Ecography (Cop) 41: 647–660. https://doi.org/10.1111/ecog.02874 MacKenzie, D. I., Nichols, J. D., Hines, J. E., Knutson, M. G. & Frank- lin, A. B. 2003: Estimating site occupancy, colonization, and local ex- tinction when a species is detected imperfectly. Ecology 84: 2200–2207 MacKenzie, D.I. & Nichols, J.D. 2004: Occupancy as a surrogate for abundance estimation. Anim Biodiversity Conservation 27.1: 461–467 Manly, B. F . L., McDonald, L., Thomas, D. L., McDonald, T. L. & Erickson, W. P . 2002: Resource Selection by Animals: Statistical Design and Analysis for Field Studies. Springer Netherlands 19/2 • 2020, 337–347 347 Kuenda Laze Identifying habitat use of Ursus arctos, Lynx lynx martinoi and Canis lupus lupus in Albanian forests using occupancy modelling May, R., Van Dijk, J., Wabakken, P ., Swenson, J. E., Linnell, J. D., Zimmermann, B., Odden, J., Pedersen, H.C., Andersen, R. & Landa, A.May R, Van Dijk J, Wabakken P , et al 2008: Habitat differentiation within the large-carnivore community of Norway’s multiple-use landscapes. Journal of Applied Ecology 45:1382–1391. https://doi. org/10.1111/j.1365-2664.2008.01527.x Melovski, D., von Arx, M., Avukatov, V., Breitenmoser-Würsten, C., Đurović, M., Elezi, R., Gimenez, O., Hoxa, B., Hristovski, S., Ivanov, G., Karamanlidis, A. A., Lanz, T., Mersini, K., Perović, A., Ramadani, A., Sanaja, B., Sanaja, P ., Schwaderer, G., Spangenberg, A., Stojanov, A., T rajçe, A. & Breitenmoser, U. 2018: Using questionnaire surveys and occupancy modelling to identify conservation priorities for the Critically Endangered Balkan lynx Lynx lynx balcanicus. Oryx 1–9. https://doi.org/10.1017/S0030605318000492 METI, NASA 2011: The Ministry of Economy, T rade, and Industry (METI) of Japan, the United States National Aeronautics and Space Administration (NASA). Advanced Spaceborne Thermal Emission and Reflection Radiometer. ASTER Global Digital Elevation Map. https:// asterweb.jpl.nasa.gov/gdem.asp. Accessed 18 Dec 2016 Molinari‐Jobin, A., Kéry, M., Marboutin, E., Marucco, F ., Zimmermann, F ., Molinari, P ., Frick, H., Fuxjäger, C., Wölfl, S., Bled, F ., Breitenmoser‐Würsten, C., Kos, I., Wölfl, M., Černe, R., Müller, O. & Breitenmoser, U. 2017: Mapping range dynamics from opportunistic data: spatiotemporal modelling of the lynx distribution in the Alps over 21 years. Animal Conservation 21: 168–180. https:// doi.org/10.1111/acv.12369 Naves, J., Wiegand, T., Revilla, E. & Delibes, M. 2003: Endangered Species Constrained by Natural and Human Factors: The Case of Brown Bears in Northern Spain. Conservation Biology 17: 1276–1289 Ordiz, A., Milleret, C., Kindberg, J., Månsson, J., Wabakken, P ., Swen- son, J. E. & Sand, H. 2015: Wolves, people, and brown bears influ- ence the expansion of the recolonizing wolf population in Scandinavia. Ecosphere 6: 1–14. https://doi.org/10.1890/ES15-00243.1 Ordiz, A., Støen, O. G., Sæbø, S., Sahlén, V., Pedersen, B. E., Kindberg, J. & Swenson, J. E. 2013: Lasting behavioural responses of brown bears to experimental encounters with humans. Journal of Applied Ecology 50: 306–314. https://doi.org/10.1111/1365- 2664.12047 Phillips, S. J., Anderson, R. P . & Schapire, R. E. 2006: Maximum entropy modeling of species geographic distributions. Ecological Modelling 190: 231–259 Piédallu, B., Quenette, P . Y., Bombillon, N., Gastineau, A., Miquel, C. & Gimenez, O. 2019: Determinants and patterns of habitat use by the brown bear Ursus arctos in the French Pyrenees revealed by occupancy modelling. Oryx 1–10. https://doi.org/10.1017/S0030605317000321 R Development Core Team 2017: R: a language and environment for statistical computing. – R Foundation for Statistical Computing, Vienna, Austria Ražen, N., Brugnoli, A., Castagna, C., Groff, C., Kaczensky, P ., Kljun, F ., Knauer, F ., Kos, I., Krofel, M., Luštrik, R., Majić, A., Rauer, G., Righetti, D. & Potočnik, H. 2016: Long-distance dispersal connects Dinaric-Balkan and Alpine grey wolf (Canis lupus) popula- tions. European journal of wildlife research 62: 137–142. https://doi. org/10.1007/s10344-015-0971-z Swenson, J. E., Wabakken, P ., Sandegren, F ., Bjärvall, A., Franzén, R. & Söderberg, A. 1995: The near extinction and recovery of brown bears in Scandinavia in relation to the bear management policies of Norway and Sweden. Wildlife Biology 1: 11–25 von Arx, M., Breitenmoser-Wuersten, C., Zimmermann, F . & Breitenmoser, U. 2004: Balkan population. In: Status and conservation of the Eurasian lynx (Lynx Lynx) in Europe in 2001 Yackulic, C. B., Chandler, R., Zipkin, E. F ., Royle, J. A., Nichols, J. D., Campbell Grant, E. H. & Veran, S. 2013: Presence- only modelling using MAXENT: when can we trust the inferences? Methods in Ecology and Evolution 4: 236–243. https://doi.org/10.1111/2041-210x.12004 Zimmermann, B., Nelson, L., Wabakken, P ., Sand, H. & Liberg, O. 2014: Behavioral responses of wolves to roads: scale- dependent ambivalence. Behavioral Ecology 25:1353–1364. https://doi.org/10.1093/beheco/aru134 Zimmermann, F ., Breitenmoser-Würsten, C. & Breitenmoser, U. 2005: Natal dispersal of Eurasian lynx (Lynx lynx) in Switzerland. Journal of Zoology 267:381–395. https://doi.org/10.1017/S0952836905007545.