Analyses of least cost paths for determining effects of habitat types on landscape permeability: wolves in Poland.

Determining ecological corridors is crucial for conservation efforts in fragmented habitats. Commonly employed least cost path (LCP) analysis relies on the underlying cost matrix. By using Ecological Niche Factor Analysis, we minimized the problems connected with subjective cost assessment or the use of presence/absence data. We used data on the wolf presence/absence in Poland to identify LCPs connecting patches of suitable wolf habitat, factors that influence patch occupancy, and compare LCPs between different genetic subpopulations. We found that a lower proportion of cities and roads surrounds the most densely populated patches. Least cost paths between areas where little dispersal takes place (i.e., leading to unpopulated patches or between different genetic subpopulations) ran through a higher proportion of roads and human settlements. They also crossed larger maximal distances over deforested areas. We propose that, apart from supplying the basis for direct conservation efforts, LCPs can be used to determine what factors might facilitate or hinder dispersal by comparing different subsets of LCPs. The methods employed can be widely applicable to gain more in-depth information on potential dispersal barriers for large carnivores.

Introduction

Landscape connectivity (i.e., the degree to which the landscape facilitates or impedes movement among resource patches) is vital to animal population survival (Taylor et al. 1993). A lack of connectivity due to habitat fragmentation may lead to modifications in the behavior of animals, for example, changes in home range size and location or changes in movement patterns (Trombulak and Frissell 2000), and Anderson and Danielson (1997) showed that corridor quality and arrangement will affect meta-population dynamics. Furthermore, Coulon et al. (2004) found evidence that gene flow in roe deer (Capreolus capreolus) was influenced by the connectivity of the landscape. Animals that have large ranges and occur in low numbers are particularly sensitive to habitat fragmentation (Noss et al. 1996). Retaining or restoring connectivity is therefore of high conservation priority (e.g., Clevenger and Waltho 2000; Kusak et al. 2009).

However, identifying the best locations for wildlife corridors is not trivial (Epps et al. 2007). Within the last decade, the use of geographic information system analyses, in particular, techniques such as “least cost path” and “friction” analyses have become increasingly popular tools (Adriaensen et al. 2003; Epps et al. 2007; Ray 2005; Schadt et al. 2002). “Costs” and “habitat suitability” are often based on informed expert opinion, leading to subjective uncertainty that might influence the final results of the model (Ray and Burgman 2006). Costs are also commonly calculated from selectivity indices derived from radio tracking or presence–absence data (e.g., Klar et al. 2008; Schadt et al. 2002). These latter data are usually based on resident animals because of the notorious scarcity of actual dispersal data. Therefore, they may reflect habitat use of resident animals, rather than costs for dispersing animals. Pre-dispersing, dispersing, and post-dispersing Iberian lynx (Lynx pardinus) show different habitat-use patterns (Palomares et al. 2000). It seems reasonable to assume that the shift in habitat use by dispersing individuals stems more from necessity than from changes in active selection (Palomares et al. 2000). Thus, using data from resident individuals might be more adequate for establishing dispersal corridors.

A further potentially biasing aspect of the use of presence/absence and to a certain extent radio-tracking data is the likelihood of “false absences” in the database (Hirzel et al. 2002b; Pearce and Boyce 2006). False absences occur when the species is present at a location but could not be detected (inaccessible area, unequal searching effort), or the habitat would be suitable but is not used for other reasons (local extinction due to historic persecution, barriers, etc.). To circumvent this problem, at least partly, Hirzel et al. (2002b) have proposed the use of a principal component analysis-based method. This method does not depend on absence data, but it still requires the presence data to be unbiased, as far as possible (Hirzel et al. 2002b). It compares the distribution of the localities where the focal species was observed to a reference set describing the whole study area. According to its use of the multidimensional space of ecological variables, they called the procedure Ecological Niche Factor Analysis (ENFA) and developed the program Biomapper (Hirzel et al. 2002a). ENFA calculates factors from eco-geographical variables (EGVs, e.g., slope, proportional area of specific habitat types, and distances to habitat features). These EGVs describe the ecological niche of a species or population. The first of the extracted factors maximizes the marginality of the species, which is defined as the ecological distance between the species optimum and the mean habitat within the reference area. Therefore, a high marginality (values close to 1) occurs if the species lives in a very particular subset of habitat type(s) relative to the reference area. The other factors describe the specialization of the species, defined as the ratio of the standard deviation of the global distribution of values of a specific ecological variable to that of the focal species. The global tolerance of the species (the inverse of the global specialization) indicates how specialized a species is, with values close to one for euryoecious (broad niche) and values close to zero for stenoecious (narrow niche) species. ENFA has been applied often with an explicit conservation aspect to a variety of species, ranging from corals to mammals (e.g., Bryan and Metaxas 2007; Oviedo and Solís 2008). The resulting habitat suitability map can, in turn, be used to determine costs (e.g., by inversing the values) and finally to conduct a least cost path analysis (Wang et al. 2008).

In the context of ecology, least cost paths (LCPs) are traditionally employed mainly to determine sites that are potentially used as dispersal routes or that should be conserved as biological corridors (e.g., Schadt et al. 2002; Epps et al. 2007). However, beyond this, comparisons between different LCPs might provide additional insights into the ecology of a species or highlight other effects of some ecological variables that are not apparent at first sight. In this paper, we propose to compare habitat features associated with different groups of LCPs as well as their source patches in order to extract further information as to which ecological variables might facilitate or hinder gene flow. This can be done even if actual dispersal data are not available. For example, a researcher might compare LCPs leading to suitable but currently uninhabited patches with LCPs leading to inhabited areas. It can be assumed that these potential corridors are not equally well suited for dispersal and ensuing gene flow, and detailed comparisons might shed light on which factors contribute to the difference. For example, the absolute or relative (i.e., per meter of path length) costs of paths may differ. While the absolute costs of a path can be expected to influence whether an individual dispersal event will be successful, the relative costs could be more relevant in a long-term perspective. Paths of high absolute but low relative costs might be too long for a single disperser; but gene flow will be possible over several generations, even if the total length of paths might be beyond the average dispersal distance of individuals. Furthermore, absolute costs are expected to correlate positively with path length. This will make comparison between paths of unequal length difficult. Wolves have been observed to disperse over more than 700 km (Mech and Boitani 2003), so that the absolute length (and thus absolute costs) of least cost paths does not seem to be a limiting factor over the scale of this study. For our study, therefore, the use of relative costs is more appropriate.

We used this approach of comparing different types of LCPs for the Polish wolf population that is particularly suitable for several reasons: Wolf censuses have been conducted regularly in Poland since 2000 (Jędrzejewski et al. 2004, 2005), so that the distribution of the species is well known. The population serves as a source for migrants into western European countries, in particular, Germany (Salvatori and Linnell 2005), but genetic data indicate that the wolf population in western Poland, in turn, heavily depends on dispersal from the north-eastern population (W. Jędrzejewski, S. Czarnomska, and co-workers, unpubl. data). Except for the eastern to north-eastern and the Carpathian population, which are in contact with the contiguous wolf range in eastern and south-eastern Europe, most regional subpopulations are now far too small to be viable over a longer period. Therefore, a good connectivity between suitable habitats is crucial for the conservation of this and other species that are heavily dependent on forest ecosystems, since suitable habitat is not necessarily the same as accessible habitat (Eigenbrod et al. 2008). Jędrzejewski et al. (2008) recently proposed a habitat suitability model for wolves in Poland based on resource selection function. It showed that there is a much more suitable habitat for wolves than currently occupied. The suitable unoccupied habitat is located in 16 patches, each covering from 400 to 25,500 km2, isolated from the present wolf range. Furthermore, genetic analyses (Pilot et al. 2006) based on mtDNA have suggested three genetic subpopulations in Poland, the differentiation of which cannot be explained by simple Euclidean distance. That study did not include samples from wolves in western Poland, but recent analyses based on microsatellite analyses suggest that these individuals are more closely related to wolves in the north-eastern part of the country than to the Carpathian or south-eastern populations (W. Jędrzejewski, S. Czarnomska, and co-workers, unpubl. data).

The main aim of our study was to evaluate how wolf colonization and gene flow can be facilitated between some patches/subpopulations or hindered between others, using a partly novel approach: comparing different sets of least cost paths that were based on a cost grid that, in turn, was based on an ENFA habitat suitability model. We hypothesize that apparently reduced gene flow between certain patches will be reflected in higher relative costs of LCPs connecting these patches. If patches were correctly classified as suitable for wolves, we would expect that they do not systematically differ in the proportion of land cover types. In contrast, LCPs associated with reduced gene flow should have a higher proportion of anthropogenically altered landscapes and a lower proportion of preferred habitat types.

Materials and methods

Brief outline

The analyses are based on a series of steps that will be explained in more detail in the following sections (see Fig. 1 for an overview). First, we conducted an Ecological Niche Factor Analysis that not only created a habitat suitability map, but, more importantly for the following analyses, also provided marginality values (see below) for all eco-geographical variables. Second, we used these values as cost values for a cost grid. Based on this cost grid, we calculated LCPs between patches of suitable wolf habitat. Third, the main analyses consisted in comparing different sets of LCPs, but also different types of wolf patches. The first two steps are methods that have to be employed so that the main analyses can be conducted. Although they need more explanation, being more technical, we would like to emphasize that the focus of this study lies in the statistical comparison of LCPs. Because we were interested in evaluating the effect of specific variables, we employed, in most cases, series of simple non-parametric tests. We decided not to conduct binary logistic regressions because these are very sensitive to collinearity of variables. Logistic regression would search for the best out of all possible models, rather than test specific hypotheses. Due to the non-normality of the data, we used non-parametric tests in most cases. Statistical analyses were conducted using R 2.10.0 (R Development Core Team 2008).

Fig. 1

Diagram of the methods employed

Habitat suitability map

We used ENFA incorporated in the program Biomapper 4.0 (Hirzel et al. 2002a, b) to calculate a habitat suitability map for wolves in the entire area of Poland. We obtained a CORINE land cover map (©EEA, Copenhagen 2000) for Poland from http://dataservice.eea.europa.eu/dataservice/. We grouped a variety of habitat types together, resulting in the following six habitat types that we converted into separate raster maps with a grid cell size of 1 km2: Arable, Forest, Meadow, Water, Wetland, and Human (including human settlements and towns). We chose the variables based on earlier studies that showed these to be of importance to wolves (Jędrzejewski et al. 2008). Preliminary analyses using a grid cell size of 250 m², as well as using coarser grids, resulted in very similar habitat suitability maps as well as in similar patterns of least cost paths (see below). We obtained data on primary (highways and international roads) and secondary roads (express roads) from the IMAGIS® Company, Warsaw. Because the presences of primary and secondary roads were highly correlated, we did not distinguish between road types in any of the further analyses (“Road”). We did not include prey density because ungulate density is correlated with forest cover, and they are common all over Poland. Also, the combined biomass of roe deer C. capreolus, red deer Cervus elaphus, and wild boar Sus scrofa never falls below 62.5 kg/km2 (calculated based on rough census maps for ungulates, T. Borowik, unpubl.; for more detailed explanations see Huck et al. 2010). We examined whether using human population density would improve the model. This variable was correlated with Human and had a less negative marginality value. We therefore did not include human population density as a separate EGV.

We used the function “circular analysis” in Biomapper to convert the raster maps for Arable, Forest, Meadow, and Human into maps representing the proportion of each habitat type in an area of 177 km2 around each grid point (still with a grid size of 1 km²). The value was constrained by the options of the program, but was closest to the 201-km2 average territory size of wolves in Białowieża forest, Eastern Poland (Jędrzejewski et al. 2007). By choosing circles corresponding to average wolf home range sizes, we ensured that the habitat suitability map (HSM) would represent suitable areas for permanent wolf populations (for a similar approach in roe deer, see, e.g., Coulon et al. 2004). However, preliminary analyses using different circle sizes gave very similar habitat suitability maps (results not shown). This indicates that the results are robust with regard to the chosen circle size. We converted the parameters Water, Wetland, and Road into distance maps (i.e., the distance to each raster point). Although ENFA also works with correlated variables, problems may occur if variables are too highly correlated. We therefore deleted Water (highly correlated to Wetland), and Road (correlated to Human), leaving the following EGVs: Arable, Human, Forest, Meadow, and Wetland. The values were Box–Cox transformed to normalize the data.

The species occurrence map consisted of wolf records collected during the National Wolf Censuses between 2000 and 2006, representing 15,670 observations (Jędrzejewski et al. 2004, 2005, for details on methods). These observations consist primarily of tracks, direct sights, howling, prey remains, and road kills, so that location errors are negligible. In comparison with other carnivores, wolves, as a group-living species, are rather conspicuous so that observation bias will be minimal. For example, even in inaccessible terrain, wolves can be detected through howling. Furthermore, all forests in Poland are divided into small forest districts and sub-compartments that are regularly checked by forestry personnel. The organizers of the census searched those areas that might have been subject to less intense monitoring by forestry staff in concentrated actions. The data are therefore likely to represent wolf occurrences without systematic bias. Preliminary data analyses using only observations that were at least 16 km apart, so that each location represented approximately one wolf pack (117 records), gave similar results. Since the reduced data set might represent variation in habitat types within home range less accurately, we used the full data set for the calculation of the final HSM. If a species lives to a large extent in suboptimal habitat, this may lead to erroneous habitat suitability calculations in ENFA, when the algorithm assumes that the median of the species' frequency distribution will represent its optimum (Braunisch et al. 2008). We therefore used the “extreme optimum adjusted median” algorithm in Biomapper 4.0 (Braunisch et al. 2008). This approach considers both the relative availability of habitat conditions, and occasions where the real optimum for the species is located at an extreme end of the environmental conditions of the study area (Braunisch et al. 2008). The number of factors included for the calculation of the HSM was recommended following the “broken stick” method (i.e., factors with eigenvalues larger than expected from randomly breaking a stick of the same total length (MacArthur 1957). The HSM was evaluated using a 20-fold cross-validation and the “continuous Boyce index” that, based on a Spearman rank correlation coefficient but using changing size classes, assesses the ability of the model to consistently predict several levels of suitability (Boyce et al. 2002; Hirzel et al. 2006). We also checked our model by comparing the percentage of wolf records on unsuitable vs. suitable or good areas within the range of the species (i.e., using a buffer with the width of an average group home-range size around all records), relative to the percentage coverage of this area.

Suitable wolf patches

In a previous study, Jędrzejewski et al. (2008) determined 33 patches over 400 km2 where wolves occurred with a probability of at least 30% (“wolf patches”). We used these wolf patches as an independently derived proxy for areas suitable for wolves. Based on the original CORINE map, we calculated for each patch the percentage area of different habitat types and of the total length of primary and secondary roads in kilometers per 100 square kilometers. Depending on whether wolves had been recorded at least three times in the National Wolf Censuses, the wolf patches were considered either as populated (patch_pop, N = 17) or unpopulated (patch_unpop, N = 16). Please note that this does not necessarily imply stable populations. Preliminary analyses suggested similar results regardless of the exact threshold for “populated”. Furthermore, we compared populated patches in eastern Poland (patch_East; patches 1–10 in Fig. 2) with populated patches in other regions (patch_other). Patches in the East are generally more densely populated by wolves. Due to the small number of southern patches, it was not feasible to use more detailed analyses for comparing other regions, for example, using only the patches in the Carpathian Mountains with the eastern patches.

Fig. 2

Suitable wolf patches within Poland and least cost paths (LCP) connecting these patches. Patches are defined following the analyses by Jędrzejewski et al. (2008), using the same number assignation for patches. Dotted lines in the south-east indicate “genetic boundaries”

For the statistical comparisons, we proceeded as follows: as background information, we compared in a Wilcoxon's signed rank test the habitat composition of wolf patches with the general habitat composition in Poland. Here, as with all other sets of similar analyses, we corrected for multiple testing using the False Discovery Rate method, which controls the average fraction of false rejections made and has a higher probability than other methods of correctly detecting real deviations between model and data (Miller et al. 2001). Using Wilcoxon's signed rank test (which is equivalent to a Mann–Whitney U test), we compared populated with unpopulated patches, and populated patches in the East with those in the rest of Poland to test the following specific hypotheses: if patches had been correctly assigned as suitable in the study by Jędrzejewski et al. (2008), we would not expect any significant differences between populated and unpopulated patches. Since wolf densities tend to be higher in eastern Poland, we expected to find for these patches either higher proportions of natural habitat (Forest and Wetland) and/or higher proportions of human modified landscapes (MaxOpen, Arable, Meadow, Human, Road).

Least cost paths

During dispersal events, animals will base moving decisions on the actual habitat rather than general suitability. We therefore used the Marginality vector of ENFA (Table 1) to assign costs to each habitat type of the CORINE map instead of converting the suitability map directly into a cost map. The marginality values were stretched between 100 (for the lowest value) and one. We then converted these values to costs by subtracting 101, using finally the absolute (rather than negative) values. Thus, Forest had a cost of one while Arable had a cost of 58. For Water, we assigned the same value as for Wetland (28), and for primary roads, we assigned the value as for Human (40), while for secondary roads we assigned half this value (Table 1). In ArcView, we combined the cost map for habitat types and for roads by simply adding the values. The combined cost map was then used to determine LCPs as described above. For each LCP, Pathmatrix (Ray 2005) provides the Euclidean distance between the patches that are connected by this path, the total length of the path, and the total cost of the path. Additionally, we divided the cost of the path by its length to obtain the “Relative Cost”. We also calculated the length (in kilometers) of path segments that crossed non-forested areas and determined the maximum value for each LCP (MaxOpen).

Table 1

Score matrix of the first four factors derived from ecological niche factor analysis, eigenvalues and the percentage of specialization explained by each factor for developing a habitat suitability map for wolves in Poland

Factor
Habitat class	1	2	3	4	cost
	(Marginality)a	(Special.1)b	(Special.2)b	(Special.3)b
Arable	−0.70	0.69	0.57	0.42	58
Human	−0.28	0.46	0.29	0.26	40
Forest	0.65	0.53	0.72	0.55	1
Meadow	0.11	0.17	0.20	0.20	24
Wetland	−0.01	0.01	0.18	0.65	28
Eigenvalue	3.25	2.68	1.46	1.09
explained specialization	35.7	29.4	16.0	11.3

The last column gives costs assigned to different habitat types according to the marginality value. For specialization factors only absolute values are shown since signs are arbitrary (Hirzel et al. 2002b)

aNegative coefficients indicate that wolves avoid areas with a high proportion of this habitat type, while positive values indicate preference.

bThe higher the value, the more restricted is the subset of values that are found in the presence of wolves

Since LCPs are one-dimensional, we created buffers of 500 m to either side of the path (i.e., a total width of 1 km). For each LCP buffer, we calculated the percentage area of the different habitat types and the total length of roads intersecting this buffer. Preliminary analyses using average patch size or the number of forest fragments per square kilometer instead of total area showed these parameters to be highly correlated. Therefore, we used only the most straightforward measure (i.e., total area) for further analyses.

As general background information, we created, for each LCP, two random routes that had the same length and the same internal angles, but different start and ending points (using the ArcView extension Alternate Routes Jenness 2004). It was not feasible to use a higher number of random paths and not necessary because of the high number of matched pairs for comparisons. We calculated the same parameters for random routes as for the original LCPs. We compared relative costs and habitat composition of LCPs with the matched average value of the random paths in a Wilcoxon's matched pair test.

Using Wilcoxon's signed rank tests (except when testing for differences in costs, see below), we compared LCPs running from populated to unpopulated wolf patches (LCP_unpop) with those running between two populated patches (LCP_pop). Paths running between two unpopulated patches or of less than 2 km length were not considered. For this analysis, we did not include LCPs “crossing genetic boundaries” (see below).

In previous studies, three genetic subpopulations of wolves have been determined in Poland, based on frequencies of mtDNA haplotypes (Pilot et al. (2006, combined with unpubl. data by W. Jędrzejewski, S. Czarnomska and co-workers) Fig. 2). We compared LCPs connecting wolf patches that were populated by wolves belonging to the same genetic subpopulation (LCP_same) with those that crossed “genetic boundaries” (LCP_cross), where wolves belonged to different genetic subpopulations. Only LCPs longer than 2 km that were connecting populated patches were considered for this analysis.

Both LCP_cross and LCP_unpop can be considered to imply reduced gene flow. We predicted that this should be reflected in higher relative costs of these paths. Using a generalized linear model (GLM) with a Gamma error distribution because of a non-normal error structure, we compared the log-transformed relative cost of LCP_cross with LCP_same, and of LCP_unpop with LCP_pop. We included the log-transformed total length of the LCPs as fixed factor because preliminary analyses correlating the relative costs of LCPs with the corresponding total length of the paths were negatively correlated (Spearman correlation using log-transformed values, R = −0.44, p << 0.001). This means that shorter LCPs were relatively more costly than long LCPs. Relative costs of random paths were not correlated with total path length (R = 0.02, p = 0.63).

Reduced gene flow might also be connected to specific landscape types. We predicted that LCP_cross and LCP_unpop should include larger proportions of anthropogenic modified landscapes (MaxOpen, Arable, Human, Road and potentially Meadow), and lower proportions of Forest and Wetland than the corresponding LCP_same and LCP_pop.

Results

Ecological niche factor analysis and habitat suitability

The four most important factors of the ENFA explained 92.4% of wolf specialization (Table 1). The score matrix of the selected factors shows that Forest (positive selection) and Arable (avoidance) are most important for explaining the marginality of the species. The Spearman correlation coefficient (“continuous Boyce index”, sensu Hirzel et al. 2006) was 0.91 ± 0.19 for a window size of 40, thus indicating good predicting power of the model. Only 21.7% of wolf records lay in areas that were considered unsuitable by ENFA, which comprised 273,481 km2 (87.6%) of Poland or 55.6% of the wolf occurrence area.

Habitat structures in wolf patches

Wolf patches differed in all parameters, but Meadow from the average conditions in Poland (Wilcoxon's matched pair test, all p < 0.001, Table 2). None of the comparisons between patch_pop and patch_unpop was significant. Only parameter Human was significant when comparing patch_East with patch_other (Wilcoxon's signed rank test, p = 0.03, not significant after correcting for multiple testing). Although it was not possible to compare the values for the Carpathian and south-eastern region in Poland (patches 11–14) with the other regions, the values given in Table 2 indicate that this area is particularly densely populated by humans.

Table 2

Median parameter values (given in percentage except for Road, which is given in kilometers per 100 square kilometers) for the entire area of Poland and for wolf patches shown for populated and unpopulated suitable wolf patches, as well as for populated patches in the East, in other regions of Poland (including the four patches in the South-East), and additionally those in the south-east separately

Parameter	Whole Poland	Patches of suitable habitat		Patches populated by wolves
		Populated	Unpopulated	East		Other	(South-East)
		N = 17	N = 16	N = 8		N = 9	N = 4
Arable	43.8a	16.0 (12.7; 20.0)	24.0 (21.1; 25.5 )	16.0 (14.4; 19.6)		19.1 (11.6; 20.0)	10.8
Forest	29.1a	55.6 (51.2; 60.8)	52.2 (50.1; 56.1)	53.0 (48.5; 62.2)		58.3 (54.9; 60.8)	47.9
Meadow	12.3	11.4 (9.6; 15.6)	10.4 (8.7; 14.7)	14.3 (10.5; 17.2)		11.4 (9.6; 11.4)	11.9
Human	8.7a	4.6 (2.6; 7.1)	6.0 (4.2; 6.9)	3.6 (1.9; 4.3)	+	7.1 (4.6; 8.1)	7.9
Wetland	4.8a	7.1 (6.0; 8.1)	6.3 (5.8; 7.7)	7.0 (5.5; 8.4)		7.1 (6.4; 7.4)	6.5
Road	17.2a	12.0 (9.7; 14.2)	13.2 (11.0; 15.1)	10.1 (7.4; 13.5)		12.5 (11.8; 18.7)	11.1

For patches, the lower and upper quantiles are given in brackets

+Significant difference (but not after correcting for multiple testing, p = 0.03) between different types of patches (to the left and right of the sign)

aSignificant difference between the coverage of that variable in Poland and the average value for all patches

As expected, LCPs (Fig. 2) had much lower relative costs than random paths of similar length. Likewise, all other parameters differed from the random value (Wilcoxon's matched pair test, all p < 0.006; Table 3).

Table 3

Parameter values (median and lower; upper quantiles, except for relative costs: mean ± 95% confidence interval) for random paths, least cost paths (LCPs) among all patches of suitable habitat for wolves (LCP_all), LCPs connecting two populated patches (LCP_pop), and LCP connecting populated patches with uninhabited patches (LCP_unpop), LCPs crossing genetic boundaries (LCP_cross), and connecting patches inhabited by the same genetic subpopulation (LCP_same)

Parameter	Random paths	p values	LCP_all	LCP_pop	p values	LCP_unpop	LCP_same	p values	LCP_cross
	N = 806		N = 403	N = 81		N = 203	N = 81		N = 56
Relative cost	39.6 (±0.6)	a	5.8 (±0.2)	5.6 (±0.2)	a	5.2 (±0.1)	5.6 (±0.2)	a	6.9 (±0.2)
MaxOpen [km]	7.6 (6.0; 9.6)	a	6.0 (4.7; 7.6)	5.0 (3.7; 6.8)		5.0 (4.6; 7.2)	5.0 (3.7; 6.8)	a	8.4 (6.0; 8.9)
Arable [%]	46.1 (40.1; 52.1)	a	12.1 (10.0; 14.1)	12.1 (9.3; 13.9)		11.8 (9.2; 13.8)	12.1 (9.3; 13.9)		12.6 (11.3; 13.9)
Forest [%]	25.3 (20.7; 30.7)	a	70.4 (66.7; 74.6)	70.5 (64.4; 73.9)	+	71.5 (67.7; 75.7)	70.5 (64.4; 73.9)		69.0 (66.1; 70.5)
Meadow [%]	11.0 (8.9; 13.2)	a	8.0 (6.5; 9.7)	8.3 (6.6; 10.4)		7.8 (6.3; 10.0)	8.3 (6.6; 10.4)		8.6 (7.9; 9.7)
Human [%]	8.0 (5.5; 11.4)	a	3.5 (1.8; 4.9)	2.3 (1.2; 3.7)	+	2.9 (1.6; 4.1)	2.3 (1.2; 3.7)	a	5.2 (4.2; 6.3)
Wetland [%]	6.0 (5.0; 7.2)	a	4.3 (3.4; 6.0)	4.6 (3.2; 7.5)		4.5 (3.4; 6.2)	4.6 (3.2; 7.5)	a	3.9 (3.4; 4.3)
Road [km/100 km path length]	8.6 (7.6; 9.8)	a	8.1 (7.0; 9.5)	7.2 (6.2; 8.6)	a	8.1 (7.0; 9.2)	7.2 (6.2; 8.6)	a	8.2 (7.3; 9.7)

P probability (Wilcoxon's matched pair and signed rank tests, see text for details)

aSignificant difference between the values to the left and right

+Significant difference (but not after correcting for multiple testing)

Taking path length into account (because relative costs decline with distance, see Methods), LCPs between two populated patches had significantly lower costs than LCPs between populated and unpopulated patches (GLM with Gamma error distribution, t = 3.17, p = 0.002, res. deviance based on 280 degrees of freedom). Likewise, LCPs between populated wolf patches that were crossing “genetic boundaries” had significantly higher relative costs than LCPs within the same genetic subpopulation (GLM with Gamma error distribution, t = 2.07, p = 0.041, residual deviance based on 133 df; Fig. 3).

Fig. 3

Correlation between (log-transformed) relative costs per meter LCP with (log-transformed) length of LCP for paths connecting same genetic subpopulations and paths crossing genetic boundaries (GLM with Gamma error distribution, t = 2.07, p = 0.041, res. deviance based on 133 df)

LCPs between two populated patches had significantly less Roads in them than LCPs leading to unpopulated patches (Wilcoxon's signed rank test, W = 10,436, p < 0.001), while, for Forest and Human the differences, were not significant after correcting for multiple testing (W = 9,536 and 10,436, p = 0.035 and 0.039 for Forest and Human, respectively, Table 3).

LCPs crossing “genetic boundaries” had significantly higher proportions of MaxOpen (W = 3450.5, p < 0.001), Human (W = 3,931, p < 0.001), Wetland (W = 1,644, p = 0.006), and Road (W = 3,146, p < 0.001) than LCPs between genetically more similar subpopulations.

Discussion

Habitat suitability and costs of paths

Our results demonstrate that detailed comparisons of LCPs and potentially suitable habitat patches, as determined by habitat suitability analysis, can offer additional insights into what kind of habitat features might hinder or facilitate dispersal of a particular species. We are not aware of any study that has compared parameters associated with different sets of LCPs. Depending on the kind of data available and the specific questions addressed in a study, some other approaches might also be appropriate. For example, Fortin et al. (2005) proposed a step selection function by which successive movement segments are compared against random segments. This approach is, in some respects, similar to our “background” comparison of LCPs with random paths. This method, however, relies on relatively continuous radio-tracking data or at least actual movement routes. In contrast, at our scale, we did not have radio-tracking data available and used complete LCPs. Furthermore, we were not primarily concerned with comparing the LCPs with random conditions—this would be rather a circular argument since the algorithm is programmed to find or avoid particular landscape types. Rather, our main concern was to find differences between different subsets of LCPs. In future studies, it might be useful to incorporate this aspect in step selection procedures. We suggest that such gained knowledge may help to inform decisions on what conservation efforts might be particularly useful. Our analysis not only provides suggestions, where biological corridors should be protected within Poland, but also points to specific problems that appear to be particularly important in hindering free dispersal of wolves and thus reducing the gene flow between populations.

The program Pathmatrix (Ray 2005) was successful in determining LCPs that were significantly less costly than random routes. Their lengths lie well within recorded dispersal distances for wolves that have been observed to cover more than 700 km (review in Mech and Boitani 2003). Nevertheless, such LCPs do not necessarily represent corridors that are actually used by animals. Proper validation of both the habitat suitability map and the LCPs will only be possible in the future when new data on colonization or dispersal events become available. In the absence of published data on real dispersal events, the model can be evaluated only indirectly; known cases of reduced dispersal (unpopulated patches, genetic differentiation) should reflect higher relative costs. Direct comparison of costs of different paths is difficult because of the negative relationship between relative costs and path length (Fig. 3). Zimmermann and Breitenmoser (2007) found a similar correlation in a corridor analysis in a Eurasian lynx population in Switzerland, without exploring that result further. The initially surprising finding can be explained by the way patches were defined a priori as suitable (i.e., including mainly habitat with low costs). Two patches that are recognized as separate will be necessarily divided by some less suitable habitat. Matrix between patches located further apart is more likely to contain some suitable habitat, so the costs of paths will not increase linearly. Paths between patches lying within a short distance from each other are unlikely to pass through much suitable habitat (otherwise, patches would not be disconnected). Future studies should keep this caveat in mind when comparing LCPs. By including path length as a factor in the GLM, we were able to show that paths that are associated with reduced gene flow (i.e., leading to uninhabited patches or connecting different genetic subpopulations of wolves) had significantly higher relative costs than other paths. With the remaining analyses, we determined what habitat types or structures might act as barriers or facilitate dispersal.

Differences among patches

When interpreting the results, it has to be kept in mind that emerging differences between the types of LCPs or patches are differences additional to those that distinguish LCPs from random paths or wolf patches from average conditions in Poland. Forest (positive selection) and Arable (negative selection) were the most important factors for explaining the marginality of the species and thus prominent in determining costs of movements. Likewise, wolf patches differed most strongly with respect to Forest and Arable from average availability of these habitats in Poland.

In confirmation to our prediction, patches in eastern Poland tended to be associated with a less dense human population and a lower density of roads than patches in other regions of the country. Neither populated nor unpopulated patches differed in any of the investigated parameters, indicating that, like in a study by Jędrzejewski et al. (2008), our model regarded all of these patches as potentially suitable for wolves. Thus, it is indeed the quality of the (potential) corridors and not the quality of the particular patch that appears to influence the connectivity to other patches. Occasionally, suitable wolf patches are effectively isolated by a very dense encirclement by settlements and towns and a dense road network. If these settlement belts are very broad, it will be difficult for animals to find gaps to enter or leave otherwise suitable patches. Thus, the particularly high human population density at the edge of the Carpathian mountains might also contribute to the genetic distinction of the wolf populations in south-eastern Poland.

Effect of habitat features on permeability

All LCPs have relatively high values for Forest and low values for arable land. Similarly, wolves studied in various populations tended to prefer forested or at least shrub habitats if available and usually avoided areas with high anthropogenic influence, even if they can adapt to these conditions to a certain degree (Ciucci et al. 2003; Jędrzejewski et al. 2004, 2005; Mech and Boitani 2003; Meriggi et al. 1991). Hence, differences are more likely to occur in those parameters that are not crucial for basic requirements. As predicted, both LCPs crossing “genetic boundaries” and LCP connecting populated with unpopulated patches of suitable wolf habitat included a higher proportion of landscapes that are strongly influenced by humans, in particular Human and Road, than LCP_same and LCP_pop. The factor Road was strongly correlated with Human, so that it is difficult to determine which of the two is more detrimental to wolf dispersal. Probably, settlements and towns are more feared by wolves, but roads might be actually more dangerous due to frequent collisions with vehicles (Lovari et al. 2007). Strong effects of urbanized areas on genetic differentiation where also found for roe deer (Wang and Schreiber 2001). The significant result for Forest (higher forest cover for LCP_unpop) does not appear to be biologically meaningful, but might rather reflect correlation-effects with other variables. A higher proportion of Wetland might contribute in two ways to facilitate dispersal. Firstly, of all “open” land cover types, wetlands are probably least disturbed by humans and may therefore be “safer” traveling routes for wild animals. Secondly, they might offer wolves additional hunting possibilities, for instance on beavers Castor fiber (W. Jędrzejewski et al., unpubl. data). In general, however, large areas of open habitat seem to hinder dispersal, as indicated by the higher maximum value of MaxOpen for LCP_cross compared with LCP_same.

Whereas it might be unfeasible to offset the barrier effect of settlements or towns (unless belts of forests and parks already exist), the effect of roads might be mitigated by wildlife passages (“green bridges” and similar structures, Jędrzejewski et al. 2009). The effectiveness of these structures for ensuring landscape connectivity not only for large carnivores, but also for ungulates and other animals has been recently demonstrated in Croatia (Kusak et al. 2009).

Conclusions

Many studies have found an influence of the environment on animal movements (Coulon et al. 2008; e.g., Fortin et al. 2005). Elk Cervus canadensis movements, for example, are influenced by various components of their environment, including the distribution of wolves (Fortin et al. 2005). Our approach differs as far as we did not directly measure movements, but deduce it from a combination of known habitat preferences (allowing the calculation of least cost paths) and genetic information (allowing distinguishing between paths with apparently more or less gene flow).

Our approach has not yet been used to compare alternative LCPs connecting roughly the same source patches. However, it could easily be used to assess which of several potential corridors merits conservation effort most or is most promising in connecting fragmented populations of different species. Furthermore, it should be suitable for a more detailed incorporation of genetic distance data.