Simulating impacts of rapid forest loss on population size, connectivity and genetic diversity of Sunda clouded leopards (Neofelis diardi) in Borneo
Population genetics Leopards Forests Heterozygosity Borneo Conservation science Species diversity Land use
Erin L. Landguth Andrew J. Hearn Yadvinder Malhi David W. Macdonald Ewan A. Macdonald Samuel A. Cushman
Habitat loss is the greatest threat to biodiversity in Borneo, and to anticipate and combat its effects it is important to predict the pattern of loss and its consequences. Borneo is a region of extremely high biodiversity from which forest is being lost faster than in any other. The little-known Sunda clouded leopard (Neofelis diardi) is the top predator in Borneo and is likely to depend critically on habitat connectivity that is currently being rapidly lost to deforestation. We modeled the effects of landscape fragmentation on population size, genetic diversity and population connectivity for the Sunda clouded leopard across the entirety of Borneo. We modelled the impacts of land use change between the years 2000, 2010 and projected forwards to 2020. We found substantial reductions across all metrics between 2000 and 2010: the proportion of landscape connected by dispersal fell by approximately 12.5% and the largest patch size declined by around 15.1%, leading to a predicted 11.4% decline in clouded leopard numbers. We also predict that these trends will accelerate greatly towards 2020, with the percentage of the landscape being connected by dispersal falling by about 57.8%, the largest patch size falling by around 62.8% and the predicted clouded leopard population falling by 62.5% between 2010 and 2020. We predicted that these large declines in clouded leopard population size and connectivity will also substantially reduce the genetic diversity of the remaining clouded leopard population.
/resource/img/favicon.ico
Figures

Abstract

Habitat loss is the greatest threat to biodiversity in Borneo, and to anticipate and combat its effects it is important to predict the pattern of loss and its consequences. Borneo is a region of extremely high biodiversity from which forest is being lost faster than in any other. The little-known Sunda clouded leopard (Neofelis diardi) is the top predator in Borneo and is likely to depend critically on habitat connectivity that is currently being rapidly lost to deforestation. We modeled the effects of landscape fragmentation on population size, genetic diversity and population connectivity for the Sunda clouded leopard across the entirety of Borneo. We modelled the impacts of land use change between the years 2000, 2010 and projected forwards to 2020. We found substantial reductions across all metrics between 2000 and 2010: the proportion of landscape connected by dispersal fell by approximately 12.5% and the largest patch size declined by around 15.1%, leading to a predicted 11.4% decline in clouded leopard numbers. We also predict that these trends will accelerate greatly towards 2020, with the percentage of the landscape being connected by dispersal falling by about 57.8%, the largest patch size falling by around 62.8% and the predicted clouded leopard population falling by 62.5% between 2010 and 2020. We predicted that these large declines in clouded leopard population size and connectivity will also substantially reduce the genetic diversity of the remaining clouded leopard population.

Citation: Macdonald EA, Cushman SA, Landguth EL, Hearn AJ, Malhi Y, Macdonald DW (2018) Simulating impacts of rapid forest loss on population size, connectivity and genetic diversity of Sunda clouded leopards (Neofelis diardi) in Borneo. PLoS ONE 13(9): e0196974. https://doi.org/10.1371/journal.pone.0196974

Editor: Luciano Bosso, Universita degli Studi di Napoli Federico II, ITALY

Received: January 10, 2018; Accepted: April 18, 2018; Published: September 12, 2018

This is an open access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.

Data Availability: All relevant data are within the paper and its Supporting Information files.

Funding: This work was funded by grants to D.W.M. by the Robertson Foundation and the Recanati-Kaplan Foundation; E.A.M. is a Kaplan Scholar and supported by the Epply Foundation and Woodspring Trust. S.A.C. was supported by the US Forest Service, Rocky Mountain Research Station during this work.

Competing interests: The authors have declared that no competing interests exist.

Introduction

To assess the impacts of habitat loss and fragmentation on populations and to develop effective conservation responses it is necessary to quantify the impacts of landscape change on population extent, density and connectivity [1, 2]. However, it is difficult to predict population connectivity reliably [2, 3] given that it is scale, species and system dependent [4]. Specifically, population connectivity is the result of the combined effects of the distribution and density of the population, composition and configuration of the landscape, and species-specific dispersal characteristics (including effects of different landscape features on movement, and the effects of sex and age differences in dispersal), and how these combine to shape the dispersal kernel [5]. In most populations there is substantial uncertainty about species distributions and densities [4], how different landscape features affect movement [6], and limited understanding of species dispersal abilities [7]. These compounding uncertainties limit the confidence that can be placed in the results of many connectivity analyses.

In addition to major gaps in knowledge about the ecology of organisms and their responses to landscape conditions, there are significant technical challenges in producing fine-scale, spatially explicit predictions of population connectivity across large areas. Traditional approaches to connectivity modelling have been based on mapping least cost paths or corridors between a few “core” habitats or protected areas [8–10], which limits the generality and scope of inferences [5]. For many questions it is better to adopt spatially synoptic modelling in which connectivity is assessed continuously across space, from all locations to all others [11–13]. The challenge arises because the tools commonly available for least cost path and corridor analysis have not enabled such synoptic analysis [14]; however, recent technical advances and improved computational power combine to enable scientists to predict more robustly the effects of landscape structure and fragmentation on population connectivity in a spatially synoptic framework. For example, resistant kernel [15–17] and factorial least-cost path approaches [18–20], coupled with landscape pattern analysis [21–23], provide a framework within which to predict the location of core habitats, barriers, and movement corridors for a range of dispersal abilities across vast landscapes [5, 24–26].

Tropical forests are among the most biodiverse habitats on earth, and provide a myriad of environmental services [27–30]. However, these habitats are being lost at a devastating rate, with estimates suggesting that approximately 230 million ha of forest were destroyed globally between 2000–2012 [31]. Of all countries, it was estimated that Indonesia had experienced the largest increase in the rate of forest loss [31], and that Borneo lost 30.2% of its forest area between 1973 and 2010 [32].

The Sunda clouded leopard (Neofelis diardi) is the apex carnivore in Borneo, an ecologically important role that adds urgency to remedying the currently minimal knowledge of its ecology, which is based largely on observations in captivity [33], anecdotal reports [34, 35] and a few recent field studies [36, 37]. However, despite this lack of ecological knowledge, the clouded leopard is charismatic and attractive to a global public [38] and has the potential to act as a flagship for wider biodiversity conservation on Borneo.

Brodie, Giordano [9] investigated clouded leopard corridors as part of a multi-species conservation approach for Sarawak. However, they employed least cost path and circuit theory to identify the optimal corridors between chosen pairs of protected areas—a completely different methodology to the spatially synoptic approaches we adopt (resistant kernel and factorial least cost path analysis; Cushman, Lewis [14]). Specifically, we address uncertainties in landscape resistance by convening a panel of experts on the species and combining their opinions. The use of expert opinion is increasing in conservation, in part due to the increasing urgency of conservation decision-making, and in particular in situations where empirically derived field data are limited [39]. We address uncertainties in dispersal ability by modelling across two levels of dispersal (S1 File), and address uncertainty in species distribution by estimating distribution and density across space as proportional to local habitat suitability as defined by the expert panel we convened. Finally, we use spatially synoptic approaches to evaluate patterns and changes in Borneo-wide connectivity for clouded leopards.

Methods

The analysis encompasses the entire island of Borneo, located in South East Asia and including territories governed by three countries (Brunei, Indonesia and Malaysia). The extent is just over 730,000km2.

We produced resistance maps for the years 2000 and 2010 for the island of Borneo based on the application of expert opinion regarding the relative habitat suitability for clouded leopards of different land cover classes. The land cover classes were taken from existing land cover maps produced for 2000 and 2010 [40] and included: water, peat swamp forest, lowland forest, lower montane forest, plantations/regrowth, lowland mosaic, montane mosaic, lowland open, montane open, urban and large scale plantation.

Expert elicitation

Experts are commonly defined as individuals who hold information about a given topic and are deferred to in its interpretation [41]. This knowledge can be acquired in a variety of ways including training, research, skills and personal experience [42]. We convened a panel of 13 experts selected on the basis of their direct involvement in field research relating to the ecology and conservation of clouded leopards.

Group approaches to expert elicitation can result in dilution of the full diversity of opinions through the dominance of one or more members of the group, polarization among subsets of members, and groupthink [41]. Panel members were therefore contacted individually and asked to provide a relative habitat suitability estimate for each land cover class. Respondents were provided with a short questionnaire that included some brief background questions regarding the duration, type and geographic location of their experience. They were then presented with a map of land uses in Borneo, along with a list of land use classes (the map, as well as the full descriptions for each class were taken from Miettinen, Shi [43], the descriptions are reproduced in Table 1). For each land use class participants were asked (i) whether they had any direct experience working in that habitat, (ii) to estimate the population density of N. diardi per 100km2 in that land use class, (iii) the relative habitat suitability of that land use class for N. diardi on an interval scale of 1–5 with 5 representing the best quality habitat and 1 representing the lowest quality habitat, (iv) for any additional information that they thought might be relevant including the sources drawn upon to produce the estimates. Respondents were asked to provide answers only for those land use classes for which they felt comfortable and were invited to contact the authors in the event that they felt any elements of the survey were unclear. Additionally a pilot version of the survey was sent to a subset of respondents for comment prior to the deployment of the final survey.

The experts showed a high level of agreement in rating habitat suitability with ICC(3,1): 0.66 and ICC(3,k): 0.96. Descriptions of land use classes from [43].

https://doi.org/10.1371/journal.pone.0196974.t001

The number of respondents who felt confident to assess habitat suitability varied between each land use class ranging from a low of 7 (Peatswamp forest, Lowland open, Large scale plantation) to a high of 12 (Lowland forest, Lower Montane forest). The frequency distributions of scores allocated to each habitat class are presented in Fig 1. The Intraclass Correlation Coefficient (ICC) was high (ICC(3,k): 0.96) and in excess of frequently quoted guidelines for interpretation of inter-rater agreement measures (>0.60 good, >0.75 excellent) [44, 45]. This provides confidence in the level of agreement between experts in assessing habitat suitability.

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

The habitat suitability score for each land use class was calculated from the mean of the habitat suitability scores provided by each respondent, this follows several similar approaches in the literature and is appropriate because responses were provided in an interval scale (Table 1)[6, 15, 46]. The mean habitat suitability scores showed a very tight relationship with the mean predicted population density for each land use class (R2 = 0.905, S1 File) providing further confidence in the habitat suitability scores.

Resistance mapping in 2000 and 2010

The habitat suitability scores were translated into relative resistances by inverting and scaling from a minimum of 1 to a maximum of 100. This translation assumes that landscape resistance is directly and inversely proportional to habitat suitability, and has a minimum value of 1 in the highest quality habitat and a maximum of 100 in the lowest quality habitat. This is a very common approach that is widely used in connectivity modelling, and habitat suitability is the most common proxy adopted for landscape resistance [6, 47]. We applied this landscape resistance classification to the land cover maps for the years 2000 and 2010 [40, 43] and then resampled to a 1000m pixel size with bilinear interpolation. The resampling was done to facilitate connectivity analysis that is extremely computationally demanding and becomes intractable at fine pixel size over large areas. Cushman and Landguth [48] found that resampling resistance maps to coarser pixel size with bilinear interpolation has relatively little effect on predicted gene flow and thus is appropriate in this case to enable connectivity modeling.

Resistance mapping for 2020

Cushman et al. [49] used a multi-scale random forest machine learning modelling approach [50] to predict deforestation risk across Borneo in 2020 from observed patterns of deforestation in relation to topography, landcover, human population density, and conservation status of lands between 2000 and 2010. This deforestation risk map provided the probability that a pixel that was forested in 2010 would be deforested by 2020. We used this deforestation risk map to translate the resistance map described above for 2010 to one predicted for 2020. This was done by multiplying the probability of deforestation by the resistance value for that location if it were plantations/regrowth and adding this net predicted change in resistance to the resistance value for that pixel in 2010. We did this because forest pixels lost to deforestation were most likely to transition either to plantations or regrowth.

Defining source points for dispersers in connectivity modelling

We used the resistance layers for 2000, 2010 and 2020 to seed locations for dispersing clouded leopards. Given that the resistance surface is an inverted and scaled version of the habitat suitability map, using the resistance surface is equivalent to seeding dispersers based on estimated habitat suitability. This was done by creating a grid of random values between 0 and 1 of the same dimension and cell size as the resistance grids and then multiplying this random grid by the inverse of the resistance value for the 2000 landscape. We chose a threshold value and selected all pixels with values larger than the threshold to use source points. We selected a threshold value that would produce 2,954 source points as this plausibly approximates a conservative estimate of the effective population size of clouded leopards on Borneo [51]; however, for the purposes of this analysis the actual number is not strictly important as we report the relative changes in connectivity and gene flow rather than the absolute values. We used this same threshold for 2010 and 2020, producing 2,595 and 2,209 source points for these years, respectively.

Connectivity modelling

Resistant kernel approach.

The resistant kernel approach to connectivity modelling is based on all-directional least-cost dispersal from a defined set of sources [15], selected as described in the section above. The resistance maps for the three dates provide resistances for all locations in the study area, in the form of the increased cost of crossing each pixel relative to the least-cost condition. These costs are used as weights in the dispersal function, such that the expected density of dispersing individuals in a pixel is down-weighted by the cumulative cost from the source, following the least-cost route [15]. The initial expected density was set to 1 in each cell containing a source point. The model calculates the expected relative density of dispersers in each pixel around the source, given the dispersal ability of the species, the nature of the dispersal function, and the resistance of the landscape [15, 16]. We used UNICOR [17] to calculate the cumulative resistant kernel density for each resistance map at 125,000 cost units (125kcu), reflecting a maximum dispersal ability of 125km in ideal habitat, and proportionately less through higher resistance landscape conditions, and at 250,000 cost units (250kcu).

Dispersal distances are not known for clouded leopards, however, a number of studies have attempted to describe allometric relationships between maximum dispersal distances and a range of physiological and life history traits for mammals in general [52, 53]. Bowman estimates that maximum dispersal distance is equal to 40*(home range size0.5) and since clouded leopards are thought to have home ranges in the region of 16km2 (N. diardi)—40km2 (N. nebulosa) we therefore estimate a likely range of dispersal distances from 160 -252km [36, 52, 54, 55]. For this study we therefore selected 125km as a plausible conservative estimate of dispersal in clouded leopards, and 250kcu as a reasonable upper limit.

Least-cost path approach.

In the second form of the connectivity modelling approach we used UNICOR [17] to produce predicted least-cost path routes from each source point to each destination point. In this factorial least cost path analysis the source points, described above, define starting and ending nodes of the least-cost paths between all pairs of individuals [3, 19]. These pair-wise least-cost paths were then combined through summation to produce maps of the density of the least-cost path network. This predicted least-cost density network was smoothed by calculating the focal mean with a 10,000m radius to better reflect local gradients of movement density [19].

Analysis of resistant kernel maps

In their raw form the resistant kernel density maps depict the expected density of dispersing individuals. It is necessary to define a threshold kernel density for delineating connected populations. We defined connected populations as those areas greater than the 10th percentile of the 2000 kernel surface and applied this threshold to all three dates (2000, 2010, 2020) to predict the extent of connected habitat across the analysis period. Given that any particular threshold is arbitrary, we decided to explore how patterns of predicted connectivity varied across various thresholds. We therefore also tested the impact of using both the 5th and 20th percentile of the 2000 kernel surface as thresholds in our analysis; the results of this are presented in S1 File. We reclassified the kernel density maps into binary form, with the value 1 where the kernel density probability of dispersal is greater than the percentile threshold and 0 where it is less.

We used FRAGSTATS [21] to calculate the percentage of the landscape, largest patch index, and number of patches that are predicted to be connected habitat for each year. The percentage of the landscape that is connected is the simplest metric of landscape composition, and quantifies how much of the study area is predicted to be connected habitat for each kernel map. The largest patch index [21] reports the extent, as a proportion of the size of the study area, of the largest patch of connected habitat. Lastly, we calculated the number of patches of internally connected habitat for each map, which is a simple measure of habitat subdivision or fragmentation.

Analysis of factorial least cost path maps

In their raw form the least cost path density maps (produced for 2000, 2010, 2020) depict the location and strength of connections uniting all pairs of source points in a least-cost network. As in the kernel analysis, it is necessary to define a threshold least-cost path density for defining meaningful connections. We chose the same threshold of the 10th percentile of the least cost path network value distribution for year 2000 and reclassified the least-cost path density maps into binary with the value 1 where the least cost path density is greater than the threshold and 0 where it is less (threshold values of 5th and 20th percentile were also tested and are presented in S1 File). We used FRAGSTATS [21] to calculate the correlation length for the least cost path density maps. Correlation length provides a measure of the average distance an organism can move within a patch before encountering the patch boundary from a random starting point [21]. The correlation length gives a global measure of the connectivity of the landscape and is a more relevant functional measure of network connectivity than more basic measures such as patch size, nearest neighbor distance, and percentage of the landscape in occupied habitat [56]. We chose to omit the other metrics since correlation length is the logical measure of the extensiveness of this network, and areal measures such as percentage of the landscape and largest patch index are less well matched to this data structure [57].

CDPOP simulations

Simulation modelling provides explicit control over pattern-process relationships [58] and enables rigorous attribution of the causes of observed patterns of population size, genetic differentiation and diversity [59–61]. By varying functional parameters, environmental characterization, and organism attributes, simulation enables scientists to investigate hypotheses about the relative influence of different factors, their interactions, and organism characteristics, such as gradients of landscape resistance, population size or dispersal ability. This provides a means for thorough evaluation of complexes of factors that would be impossible to investigate directly in the field. In addition, simulations are particularly important in evaluating scenarios of future landscape change, such as predicted patterns of deforestation, on landscape resistance and its effects on future population genetic structure.

We used CDPOP 1.0 [20] to simulate the processes of mating and dispersal as functions of the spatial patterns of resistance across the three resistance maps corresponding to the 2000, 2010 and 2020 dates at the dispersal distances of 125kcu and 250 kcu. CDPOP is an individual based, spatially explicit, landscape genetic program that simulates birth, death, mating and dispersal of individuals in complex landscapes as probabilistic functions of movement cost among them. In each of the three landscape resistance maps, we used the same source points described above in the resistant kernel analysis as locations of individual clouded leopards. We stipulated the population to have 30 loci, with 10 alleles per locus, initially randomly assigned among individuals and a mutation rate of 0.0005, which are standard parameters widely used in landscape genetics simulation modelling We used an inverse square mating and dispersal probability function, with maximum dispersal cost-weighted distance of 125kcu, matching the dispersal ability in the resistant kernel analysis. The number of offspring was based on a Poisson probability draw with mean of 2. For each of the six resistance scenarios, we ran 10 Monte Carlo runs in CDPOP to assess stochastic variability. We simulated gene flow among these locations for 200 non-overlapping generations. Past studies have shown that this is sufficient time to ensure spatial genetic equilibrium [62]. We extracted several global measures of population genetic structure for the full study area at generation 200. These include total number of alleles in the population, observed and expected heterozygosity. We analyzed the differences in these global measures of genetic structure between 2000, 2010, and 2020.

sGD prediction of local genetic diversity

We utilized the sGD program [59] to map local total number of alleles and observed heterozygosity. sGD uses a variable radius focal window to calculate locally centered values of a number of population genetic parameters, such as allelic richness, expected and observed heterozygosity, across spatially structured populations. We calculated the local number of alleles, observed heterozygosity and expected heterozygosity within local windows with width equal to 1/2 of the range of significant spatial genetic autocorrelation of the Mantel correlogram across the average of the 10 runs of CDPOP for each scenario [23, 63].

Several of the approaches used in this study are computationally demanding (e.g the spatially synoptic analysis of least cost paths, and the individual based population genetics models) and may not be practical for many conservation practitioners. We therefore tested whether it was possible to take a pragmatic shortcut in estimating landscape scale changes in allelic richness and heterozygosity by simply using the maps of landscape resistance, path density and kernel density. We used linear regression (using the lm library in R) to construct linear models relating these measures of genetic structure to cumulative resistant kernel value, factorial least cost path density and local average landscape resistance within a focal window of 10,000m radius to explore the relationship between these genetic attributes and landscape patterns of population connectivity.

Results

Landscape resistance change between 2000 and 2010 was driven by loss of forest cover along the edges of forested patches in Kalimantan and by a more diffuse pattern of perforation and attrition of forest patches in Sarawak and Sabah (Fig 2a and 2b). The change in landscape resistance predicted between 2010 and 2020 extended and accelerated the pattern observed between 2000 and 2010. Specifically, in Kalimantan additional forest loss was predicted primarily along the edges of existing forest patches while in Malaysian Borneo we predicted large changes in regional resistance as a result of predicted extensive and diffuse forest loss across broad areas (Fig 2c).

Landscape resistance maps are scaled from 1 to 100, cumulative resistant kernel maps and factorial least cost path density maps are scaled from low density in blue to high density in red.

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

Analysis of resistant kernel maps

Fig 2d, 2e and 2f shows cumulative resistant kernel density of predicted frequency of dispersal given a 125kcu dispersal threshold for the years 2000, 2010, and 2020. Change in resistant kernel density at the 250kcu threshold is shown in Figure B in S1 File.

The extent of the landscape connected by dispersal changed markedly between 2000 and 2020 (Table 2). The percentage of connected landscape declined by 12.5% between 2000 and 2010; however, this loss accelerated and there was a 57.83% decrease in the extent of the landscape connected by dispersal between 2010 and 2020 (Table 2). The relative effect of the different percentile thresholds and dispersal abilities was much less than that due to landscape change between dates (S1 File).

https://doi.org/10.1371/journal.pone.0196974.t002

Largest patch index followed a similar pattern with the extent of the largest patch of connected habitat as percentage of the landscape falling by 15.1% between 2000 and 2010, and then 62.75% between 2010 and 2020 (Table 2). As with percentage of the landscape connected by dispersal, the change in largest patch index between years had a much greater effect than choice of dispersal threshold or percentile threshold on largest patch index (S1 File).

The number of isolated patches of internally connected habitat increased linearly, from 2 in 2000, to 4 in 2010 and to 6 in 2020 (Table 2) suggesting that large connected patches were broken up and fragmented. Selecting different kernel percentile thresholds, and dispersal distances affected the trend in the number of isolated patches, highlighting the complexity of the impacts of forest fragmentation (S1 File).

Analysis of factorial least cost path maps

Fig 2g, 2h and 2i shows least cost path density for the years 2000, 2010, and 2020. The correlation length of the connected least cost path network decreased by 2.5% between 2000 and 2010 and by 9.7% between (Table 2). Correlation length was also sensitive to choice of percentile threshold, declining 9.2% as the threshold was moved from 5th to 10th percentile, and another 4.6% when moved from 10th to 20th percentile, for the year 2000 connectivity map (S1 File).

CDPOP simulation results

Over the three CDPOP scenarios (2000, 2010, 2020) there was considerable variation in simulated population size, expected heterozygosity, observed heterozygosity, average mating distance and average dispersal distance (Fig 3, Table 3). There was relatively little difference between the 125kcu and 250kcu dispersal scenarios in these population characteristics (Fig 3, Table 3). In contrast, there was a large difference between the years for each of these parameters. For example, the simulations predicted that after 200 simulated generations the population supported by the 2010 landscape was 11.4% and 14.4% smaller than that supported by the landscape in 2000, for the 125kcu and 250kcu scenarios, respectively. The simulations predicted that the clouded leopard populations in the 2020 landscape would be 62.6% and 49.8% smaller than in 2010, for the 125kcu and 250kcu scenarios, respectively.

All values represent the equilibrium position after a simulated 200 generations of mating under each landscape condition. Each simulation was run 10 times and the points indicate the mean result and the bars indicate the standard deviation.

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

https://doi.org/10.1371/journal.pone.0196974.t003

After 200 generations, there was relatively little predicted change in the equilibrium values for observed heterozygosity and inbreeding coefficient between the 2000 and 2010 landscapes, at either dispersal threshold. However, the simulations predicted that the 2020 landscape would lead to a 48.5% and 36.8% reduction in observed heterozygosity compared to 2010, for the 125kcu and 250kcu scenarios, respectively (Fig 3, Table 3). This corresponded to increases in simulated inbreeding coefficient of 130% and 225% respectively (Fig 3, Table 3).

Relationships between local genetic diversity and landscape connectivity

We explored the relationships between local average number of alleles per locus (Fig 4) and local heterozygosity (Fig 4) and focal mean landscape resistance, cumulative density of least cost paths, and cumulative resistant kernel density. Visual inspection of the scatterplots suggested little relationship between allelic richness or heterozygosity and average landscape resistance within a 10km radius focal window. Cumulative density of least cost paths and cumulative resistant kernel value evinced curvilinear relationships with both allelic richness and heterozygosity in local genetic neighborhoods.

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

Our regression analysis indicated that cumulative resistant kernel density was the dominant predictor, with an inverse transformation to linearize the relationship for number of alleles and the heterozygosity at 125kcu dispersal distance. We evaluated all combinations of local resistance, cumulative least cost path density and cumulative resistant kernel value for both the average alleles per locus and observed heterozygosity. We selected the single variable model with cumulative resistant kernel value as the only predictor variable, based on the observation that the single variable resistant kernel model was within 1% of the total deviance explained by the global model, and that none of the other variables’ single variable models explained more than about 15% of the deviance. This indicates that neither local landscape resistance nor cumulative least cost path density are good predictors of genetic diversity in a population. On the other hand, the resistant kernel models were remarkably good, explaining over 52% of deviance for both alleles and heterozygosity in the 125kcu dispersal scenario (Fig 5; this rose to excess of 80% of the deviance for both variables in the 250kcu dispersal distance scenario (S1 File)).

The fitted regression equation is shown in red overlay, and the equation and deviance explained are displayed below each scatterplot.

https://doi.org/10.1371/journal.pone.0196974.g005

Discussion

Ongoing rapid landscape change on Borneo is reducing the extent and increasing the fragmentation of habitat for forest dependent species, such as the clouded leopard. This reduction in habitat amount and connectivity has important repercussions for both genetic diversity and projected population sizes. Our analysis suggested rapidly accelerating reductions in the extent and connectivity of forest habitats, with large predicted decreases in clouded leopard population size, accompanied by increasing inbreeding coefficient and decreasing genetic diversity. This suggests that shrinkage and fragmentation of habitat across Borneo will not only result in smaller clouded leopard populations, but also populations that are potentially more vulnerable to demographic stochasticity and inbreeding depression as well as being less able to adapt to future environmental change [64]. These trends were robust across dispersal distances (125km and 250km) and across different percentile thresholds for clouded leopard core areas.

Relationships between genetic diversity and landscape connectivity

This is the first study to model effects of ongoing deforestation on population decline and changes in population connectivity and genetic diversity of any species at this scale. The spatially synoptic analysis of least cost paths across Borneo in particular was extremely computationally demanding, requiring the calculation of all least cost paths between all source points (between 715,806 and 4,114,146 least cost paths for each scenario). In addition, parameterizing and running an individual based population genetics model such as CDPOP across such a vast scale and large population size represents a considerable challenge for most practitioners. Therefore, we tested whether it was possible to take a pragmatic shortcut in estimating landscape scale changes in allelic richness and heterozygosity by simply using the maps of landscape resistance, path density and kernel density that we had already computed. Our regression analysis showed that the full model (using all three predictor variables) allowed us to explain about 60% of the deviance for both allelic richness and heterozygosity. However in all cases, the single variable model (including only kernel density) was the most parsimonious solution and explained an almost identical proportion of the deviance. Future studies that are limited by time or computational power might therefore be able to rely on cumulative resistant kernel density as a proxy for changes in measures of genetic diversity.

Conservation implications and future research

Despite a recent increase in research effort, the Sunda clouded leopard on Borneo remains elusive and one of the world’s least known felids, and this hampers the development of coherent conservation actions [51]. Nonetheless, this study highlights the severe impacts of ongoing fragmentation and loss of forest habitats on Borneo. Our analyses predicted a marked decline in both the overall area of clouded leopard core areas and a marked reduction in patch sizes with several patches being lost all together, coupled with dramatic declines in allelic richness and heterozygosity. Future conservation effort should focus on slowing the rate of land use change in Borneo, while prioritizing conservation in areas that maintain the extent of core populations and the connectivity between them. Additionally, while a small number of studies have successfully made local estimates of clouded leopard density [36, 37, 65], future research should focus on better estimates of clouded leopard abundance and density across Borneo as well as improving knowledge of the species’ basic ecology (particularly its dispersal) and how it responds to anthropogenic changes in its environment. This would improve our ability to reliably map both habitat suitability and landscape resistance. We are currently in the process of refining the models presented in this paper using empirical field data from the relatively well-studied region of Sabah, however Borneo wide initiatives are also required.

Further implications

The implications of this study extend beyond the conservation of the Sunda clouded leopard. Firstly, while our models were parameterized specifically for clouded leopards based on expert opinion they are likely to have relevance to any other medium sized forest dependent mammals on Borneo (such as the sun bear (Helarctos malayanus)) that may have similar dispersal distances and habitat associations.

Macdonald, Burnham [66] have highlighted the fact that many species (for example felids and primates) face the same threats in the same locations and may benefit from many of the same solutions. Building on this notion, [67] and [38] have shown that the clouded leopard is not only regarded as charismatic by the general public but that it has strong potential to act as an Ambassador for conservation marketing campaigns. This supports the widespread approach of protecting patches or corridors that support one or several wide-ranging, often large-bodied, species on the assumption that conservation of these umbrella species will also facilitate conservation of smaller or less mobile organisms (e.g. [68, 69] but see [5, 25]). As a generalization, Beier et al. [70] argued that habitat corridors may more effectively protect regional biodiversity if they are developed to support the movement of multiple species simultaneously, because so many species are threatened by fragmentation, a proposition which underlies the analyses presented by [9] for a basket of SE Asian species including clouded leopards.

Constraints and limitations of the approach

Some parts of our model are populated by data derived from expert opinion, and although our team includes researchers with many years of experience conducting field research on clouded leopards, maps of habitat suitability derived from expert opinion have been shown to be less reliable than empirically derived data [71–73]. In this study, the ICC suggests that survey respondents showed a high level of agreement, however inevitably the results for some land uses were more certain than others. Frequency distributions for the survey responses are shown in Fig 1, and these show that while reviewers were reasonably unanimous in their view that lowland forests represented high quality habitat, they were less consistent in their view towards peat swamp forest. The comparative uncertainty in these areas could influence the reliability of our model and highlights the need for greater research in these less understood habitat types. Reassuringly however, the estimated population density provided by the survey respondents (1.12 ind. per 100km2) falls within the empirically estimated range of population density for clouded leopards in the peat swamp forest of Sabangau (0.72 to 4.41 ind. per 100km2)[65]. For this reason it will be important to collect empirical data on population distribution, abundance, genetic diversity and gene flow to verify and improve our predictions.

Additionally, our method for estimating the landscape resistance in 2020 has inherent uncertainty. Our method involves multiplying pixel values by the probability of deforestation in each cell, and thus if we have two pixels that both have a 50% probability of deforestation by 2020, then both cells will be multiplied by 50% whereas in reality we might expect that one cell might be completely deforested and the other might remain intact. One potential solution to this inaccuracy would have been to allocate a number of cells for deforestation based on their probability of loss. Under this approach, however, we would still have no objective way of knowing which pixel would be deforested. Our approach cannot therefore be used to interpret fine scale impacts of local land use changes, although the impacts at a wider landscape scale remain robust. Nonetheless, the present urgency of guiding conservation in Borneo requires prompt action to address an unfolding biodiversity crisis. Therefore, we believe the governments of Borneo should work with conservationists, scientists and economists to utilize the best information currently available to develop proactive strategies to slow the dramatic pace of forest loss and biodiversity impacts. It was our goal in this paper to help provide some information for this effort.

Conclusion

Deforestation is on-going and accelerating across Borneo with immense implications for biodiversity. Our analysis is the first to use robust predictions of future landscape change at such a broad scale to simulate changes in population size, connectivity, gene flow and genetic diversity. Our results suggest large recent declines in clouded leopard population size and genetic diversity and predict much larger imminent declines as deforestation continues across Borneo, particularly in Sabah and Sarawak, where the pattern of forest loss is much more damaging due to its widespread, diffuse nature. Our analyses provide plausible inferences about how recent and future landscape change may affect clouded leopard populations (and probably those of other species with similar biology) across Borneo. Additional work using extensive occurrence, genetic and movement data is needed to verify and improve our predictions. However, given the immediate and urgent need for effective conservation action in the part of the world experiencing the most rapid loss of tropical forest, we hope our analyses could be of immediate utility in identifying (1) population core areas for priority conservation, (2) linkages between these core areas, (3) relationships between habitat extent, landscape resistance, population size, dispersal ability and genetic diversity to help managers and planners anticipate the effects of a range of alternative conservation scenarios. The tools we have presented here could be of use in developing and evaluating the specific impacts of alternative scenarios on population size, connectivity and genetic diversity.

Supporting information

S1 File. Supporting information for simulating impacts of rapid forest loss on population size, connectivity and genetic diversity of Sunda clouded leopards (Neofelis diardi) in Borneo.

Fig A. Plot of the mean respondent scores for habitat suitability vs the mean estimated population density for clouded leopards. The high of R2 = 0.905 shows high consistency in the estimates provided by the panel of experts. Fig B. Comparison of resistant kernel maps for two potential dispersal distances. 125kcu (as shown in Fig 2 -,A, B, C, and 250kcu, D, E, F, for years 2000 (A and D), 2010 (B and E) and 2020 (C and F). Fig C. Change in the percentage of the landscape connected by dispersal for 2000, 2010, and 2020 across varying thresholds for selection and dispersal distance. Fig D. Change in the extent of the largest patch of connected habitat as percentage of the landscape across the 5th, 10th and 20th percentiles of the cumulative kernel surface in year 2000 for 2000, 2010, and 2020 at (a) 125kcu and (b) 250kcu dispersal thresholds. Fig E. Change in the number of isolated patches of connected habitat across the 5th, 10th and 20th percentiles of the cumulative kernel surface in year 2000 for years 2000, 2010, and 2020 at (a) 125kcu and (b) 250kcu dispersal thresholds. Fig F. Changes in the correlation length of the factorial least cost path network across three density thresholds (5th percentile, 10th percentile and 20th percentile of the least cost path density in 2000), and across the three dates (2000, 2010, 2020). Fig G. Scatterplots and fitted LOWESS splines for relationship between average number of alleles per locus and focal mean landscape resistance within a 10km radius (column 1), cumulative density of least cost paths (column 2), and cumulative resistant kernel density (column 3), across the two dispersal distance scenarios (125kcu, row 1; 250kcu, row 2). Fig H. Scatterplots and fitted LOWESS splines for relationship between observed heterozygosity and focal mean landscape resistance within a 10km radius (column 1), cumulative density of least cost paths (column 2), and cumulative resistant kernel density (column 3), across the two dispersal distance scenarios (125kcu, row 1; 250kcu, row 2). Fig I. Scatterplots of average number of alleles per locus in a local neighborhood relative to resistant kernel value for a) the 125kcu scenario and b) the 250kcu scenario, and heterozygosity of the local neighborhood relative to resistant kernel value for c) the 125kcu scenario and d) the 250kcu scenario. The fitted regression equation is shown in red overlay, and the equation and deviance explained are displayed below each scatterplot.

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

(DOCX)

S2 File. Core Areas and Corridors data for PLOS—LCP and Resistance.

This .zip file contains the following shapefiles:

Landscape Resistance map for the year 2000 (“Resistance_2000.tif”)

Landscape Resistance map for the year 2010 (“Resistance_2010.tif”)

Landscape Resistance map for the year 2020 (“Resistance_2020.tif”)

Least Cost Path map for the year 2000 (“LCP_2000.tif”)

Least Cost Path map for the year 2000 (“LCP_2010.tif”)

Least Cost Path map for the year 2000 (“LCP_2020.tif”).

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

(ZIP)

S3 File. Core Areas and Corridors data for PLOS—Kernels.

This .zip file contains the following shapefiles:

Resistant kernel map for year 2000 with dispersal threshold of 125kcu (“kernel_2000_125kcu.tif”)

Resistant kernel map for year 2010 with dispersal threshold of 125kcu (“kernel_2010_125kcu.tif”)

Resistant kernel map for year 2020 with dispersal threshold of 125kcu (“kernel_2020_125kcu.tif”)

Resistant kernel map for year 2000 with dispersal threshold of 250kcu (“kernel_2000_250kcu.tif”)

Resistant kernel map for year 2010 with dispersal threshold of 250kcu (“kernel_2010_250kcu.tif”)

Resistant kernel map for year 2020 with dispersal threshold of 250kcu (“kernel_2020_250kcu.tif”).

https://doi.org/10.1371/journal.pone.0196974.s003

(ZIP)

Acknowledgments

The authors would like to thank everybody who has contributed to this manuscript. This work was funded by grants to D.W.M. by the Robertson Foundation and the Recanati-Kaplan Foundation; E.A.M. undertook this work while a Kaplan Scholar and supported by the Epply Foundation and Woodspring Trust. S.A.C. was supported by the US Forest Service, Rocky Mountain Research Station during this work. The authors would also like to thank Professors E.J. Milner-Gulland and L. Boitani who, as EAM’s thesis examiners, and along with anonymous reviewers, also contributed greatly.

References

1. Bennett AF. Linkages in the landscape: the role of corridors and connectivity in wildlife conservation: Iucn; 1999.

2. Rudnick DA, Ryan SJ, Beier P, Cushman SA, Dieffenbach F, Epps CW, et al. The role of landscape connectivity in planning and implementing conservation and restoration priorities. Issues in Ecology. 2012;16.

3. Cushman SA, McRae B, Adriaensen F, Beier P, Shirley M, Zeller KA. Biological corridors and connectivity. In: Macdonald DW, Willis KJ, editors. Key Topics in Conservation Biology 22013.

4. Cushman SA. Effects of habitat loss and fragmentation on amphibians: A review and prospectus. Biological Conservation. 2006;128(2):231–40. http://dx.doi.org/10.1016/j.biocon.2005.09.031.

5. Cushman SA, Landguth EL. Multi-taxa population connectivity in the Northern Rocky Mountains. Ecological Modelling. 2012;231:101–12.

6. Zeller KA, McGarigal K, Whiteley AR. Estimating landscape resistance to movement: a review. Landscape Ecol. 2012;27(6):777–97.

7. Elliot NB, Cushman SA, Loveridge AJ, Mtare G, Macdonald DW. Movements vary according to dispersal stage, group size and rainfall: the case of the African lion. Ecology. 2014.

8. Beier P, Spencer W, Baldwin RF, McRae BH. Toward Best Practices for Developing Regional Connectivity Maps. Conservation Biology. 2011;25(5):879–92. pmid:21797924

9. Brodie JF, Giordano AJ, Dickson B, Hebblewhite M, Bernard H, Mohd-Azlan J, et al. Evaluating Multispecies Landscape Connectivity in a Threatened Tropical Mammal Community. Conservation Biology. 2014:n/a–n/a. pmid:25065425

10. Kabir M, Hameed S, Ali H, Bosso L, Din JU, Bischof R, et al. Habitat suitability and movement corridors of grey wolf (Canis lupus) in Northern Pakistan. PloS one. 2017;12(11):e0187027. pmid:29121089

11. Cushman SA, Lewis JS, Landguth EL. Evaluating the intersection of a regional wildlife connectivity network with highways. 2013.

12. Khosravi R, Hemami MR, Cushman SA, Thuiller W. Multispecies assessment of core areas and connectivity of desert carnivores in central Iran. Diversity and Distributions. 2018;24(2):193–207.

13. Reed GC, Litvaitis JA, Callahan C, Carroll RP, Litvaitis MK, Broman DJA. Modeling landscape connectivity for bobcats using expert-opinion and empirically derived models: how well do they work? Animal Conservation. 2017;20(4):308–20.

14. Cushman SA, Lewis JS, Landguth EL. Why Did the Bear Cross the Road? Comparing the Performance of Multiple Resistance Surfaces and Connectivity Modeling Methods. Diversity. 2014;6(4):844–54.

15. Compton BW, McGarigal K, Cushman SA, Gamble LR. A Resistant-Kernel Model of Connectivity for Amphibians that Breed in Vernal Pools. Conservation Biology. 2007;21(3):788–99. pmid:17531056

16. Cushman SA, Chase M, Griffin C. Mapping Landscape Resistance to Identify Corridors and Barriers for Elephant Movement in Southern Africa. In: Cushman S, Huettmann F, editors. Spatial Complexity, Informatics, and Wildlife Conservation: Springer Japan; 2010. p. 349–67.

17. Landguth EL, Hand BK, Glassy J, Cushman SA, Sawaya MA. UNICOR: a species connectivity and corridor network simulator. Ecography. 2012;35(1):9–14.

18. Cushman SA, Compton BW, McGarigal K. Habitat Fragmentation Effects Depend on Complex Interactions Between Population Size and Dispersal Ability: Modeling Influences of Roads, Agriculture and Residential Development Across a Range of Life-History Characteristics. In: Cushman S, Huettmann F, editors. Spatial Complexity, Informatics, and Wildlife Conservation: Springer Japan; 2010. p. 369–85.

19. Cushman SA, McKelvey KS, Schwartz MK. Use of empirically derived source-destination models to map regional conservation corridors. Conservation biology: the journal of the Society for Conservation Biology. 2009;23(2):368–76. Epub 2008/11/20. pmid:19016821.

20. Landguth EL, Cushman SA. cdpop: A spatially explicit cost distance population genetics program. Molecular Ecology Resources. 2010;10(1):156–61. pmid:21565001

21. McGarigal K, Cushman SA, Ene E. FRAGSTATS v4: Spatial Pattern Analysis Program for Categorical and Continuous Maps. Computer software program produced by the authors at the University of Massachusetts, Amherst. http://www.umass.edu/landeco/research/fragstats/fragstats.html. 2012.

22. Cushman SA, Landguth EL, Flather CH. Evaluating the sufficiency of protected lands for maintaining wildlife population connectivity in the U.S. northern Rocky Mountains. Diversity and Distributions. 2012;18(9):873–84.

23. Wasserman TN, Cushman SA, Shirk AS, Landguth EL, Littell JS. Simulating the effects of climate change on population connectivity of American marten (Martes americana) in the northern Rocky Mountains, USA. Landscape Ecol. 2012;27(2):211–25.

24. Cushman S, Landguth EL. Ecological associations, dispersal ability, and landscape connectivity in the northern Rocky Mountains. 2012.

25. Cushman SA, Landguth EL, Flather CH. Evaluating population connectivity for species of conservation concern in the American Great Plains. Biodiversity and Conservation. 2013;22(11):2583–605.

26. Hand B, Cushman S, Landguth E, Lucotch J. Assessing multi-taxa sensitivity to the human footprint, habitat fragmentation and loss by exploring alternative scenarios of dispersal ability and population size: a simulation approach. Biodiversity and Conservation. 2014;23(11):2761–79.

27. Pimm SL, Raven P. Biodiversity: Extinction by numbers. Nature. 2000;403:843–5. pmid:10706267

28. Pielke RA, Marland G, Betts RA, Chase TN, Eastman JL, Niles JO, et al. The influence of land-use change and landscape dynamics on the climate system: relevance to climate-change policy beyond the radiative effect of greenhouse gases. Philosophical Transactions of the Royal Society of London Series A: Mathematical, Physical and Engineering Sciences. 2002;360(1797):1705–19. pmid:12460493

29. van der Werf GR, Morton DC, DeFries RS, Olivier JGJ, Kasibhatla PS, Jackson RB, et al. C02 emissions from forest loss. Nature Geoscience. 2009;2.

30. Foley JA, Defries R, Asner GP, Barford C, Bonan G, Carpenter SR, et al. Global consequences of land use. Science. 2005;309(5734):570–4. Epub 2005/07/26. pmid:16040698.

31. Hansen MC, Potapov PV, Moore R, Hancher M, Turubanova SA, Tyukavina A, et al. High-resolution global maps of 21st-century forest cover change. Science. 2013;342(6160):850–3. Epub 2013/11/16. pmid:24233722.

32. Gaveau DLA, Sloan S, Molidena E, Yaen H, Sheil D, Abram NK, et al. Four Decades of Forest Persistence, Clearance and Logging on Borneo. PloS one. 2014;9(7):e101654. pmid:25029192

33. Selous EM, Banks E. The clouded leopard in Sarawak. J Sarawak Museum. 1935;9:263–6.

34. Santiapillai C, Ashby KR. THE CLOUDED LEOPARD IN SUMATRA. Oryx. 1988;22(1):44–5.

35. Rabinowitz A, Andau P, Chai PP. The clouded leopard in Malaysian Borneo. Oryx. 1987;21(02):107–11.

36. Hearn AJ, Ross J, Pamin D, Bernard H, Hunter L, Macdonald DW. Insights into the the spatial and teporal ecology of the sunda clouded leopard Neofelis Diardi. Raffles Bulletin of Zoology. 2013;61(2):871–5.

37. Brodie J, Giordano AJ. Density of the Vulnerable Sunda clouded leopard Neofelis diardi in a protected area in Sabah, Malaysian Borneo. Oryx. 2012;46(3):427–30.

38. Macdonald EA, Burnham D, Hinks AE, Dickman AJ, Malhi Y, Macdonald DW. Conservation inequality and the charismatic cat: Felis felicis. Global Ecology and Conservation. 2015;3(0):851–66. http://dx.doi.org/10.1016/j.gecco.2015.04.006.

39. Kuhnert PM. Four case studies in using expert opinion to inform priors. Environmetrics (special issue). 2011;22:662–74.

40. Miettinen J, Shi C, Tan WJ, Liew SC. 2010 land cover map of insular Southeast Asia in 250-m spatial resolution. Remote Sensing Letters. 2012;3(1):11–20.

41. Martin TG, Burgman MA, Fidler F, Kuhnert PM, Low-Choy S, McBride M, et al. Eliciting Expert Knowledge in Conservation Science. Conservation Biology. 2012;26(1):29–38. pmid:22280323

42. Burgman M, Carr A, Godden L, Gregory R, McBride M, Flander L, et al. Redefining expertise and improving ecological judgment. Conservation Letters. 2011;4(2):81–7.

43. Miettinen J, Shi C, Liew SC. Deforestation rates in insular Southeast Asia between 2000 and 2010. Global Change Biology. 2011;17(7):2261–70.

44. Cicchetti DV. Guidelines, criteria, and rules of thumb for evaluating normed and standardized assessment instruments in psychology. Psychological Assessment. 1994;6(4):284–90.

45. Shrout PE, Fleiss JL. Intraclass correlations: uses in assessing rater reliability. Psychological bulletin. 1979;86(2):420. pmid:18839484

46. Siegel S, Castellan NJ Jr.. Nonparametric Statistics for the Behavioural Sciences: McGraw-Hill Book Company; 1988.

47. Mateo-Sánchez M, Balkenhol N, Cushman S, Pérez T, Domínguez A, Saura S. A comparative framework to infer landscape effects on population genetic structure: are habitat suitability models effective in explaining gene flow? Landscape Ecol. 2015;30(8):1405–20.

48. Cushman SA, Landguth EL. Scale dependent inference in landscape genetics. Landscape Ecol. 2010;25(6):967–79.

49. Cushman SA, Macdonald EA, Landguth EL, Malhi Y, Macdonald DW. Multiple-scale prediction of forest loss risk across Borneo. Landscape Ecol. in press.

50. Evans JS, Cushman SA. Gradient modeling of conifer species using random forests. Landscape Ecol. 2009;24(5):673–83.

51. Hearn AJ, Ross J, Brodie J, Cheyne SM, Haidir IA, Loken B, et al. Neofelis diardi. The IUCN Red List of Threatened Species 2015 2015.

52. Bowman J, Jaeger JAG, Fahrig L. DISPERSAL DISTANCE OF MAMMALS IS PROPORTIONAL TO HOME RANGE SIZE. Ecology. 2002;83(7):2049–55.

53. Whitmee S, Orme CDL. Predicting dispersal distance in mammals: a trait-based approach. Journal of Animal Ecology. 2013;82(1):211–21. pmid:22924343

54. Grassman LI, Tewes ME, Silvy NJ, Kreetiyutanont K. ECOLOGY OF THREE SYMPATRIC FELIDS IN A MIXED EVERGREEN FOREST IN NORTH-CENTRAL THAILAND. Journal of Mammalogy. 2005;86(1):29–38.

55. Austin SC, Tewes ME, Grassman LI Jr., Silvy NJ. Ecology and conservation of the leopard cat Prionailurus bengalensis and clouded leopard Neofelis nebulosa in Khao Yai National Park, Thailand. Acta Zoologica Sinica. 2007;53(1):1–14.

56. Cushman SA, Shirk A, Landguth EL. Separating the effects of habitat area, fragmentation and matrix resistance on genetic differentiation in complex landscapes. Landscape Ecol. 2012;27(3):369–80.

57. Short Bull RA, Cushman SA, Mace R, Chilton T, Kendall KC, Landguth EL, et al. Why replication is important in landscape genetics: American black bear in the Rocky Mountains. Molecular Ecology. 2011;20(6):1092–107. pmid:21261764

58. Epperson BK, McRae BH, Scribner KIM, Cushman SA, Rosenberg MS, Fortin M-J, et al. Utility of computer simulations in landscape genetics. Molecular Ecology. 2010;19(17):3549–64. pmid:20618894

59. Shirk AJ, Cushman SA. sGD: software for estimating spatially explicit indices of genetic diversity. Molecular Ecology Resources. 2011;11(5):922–34. pmid:21679313

60. Wasserman TN, Cushman SA, Schwartz MK, Wallin DO. Spatial scaling and multi-model inference in landscape genetics: Martes americana in northern Idaho. Landscape Ecol. 2010;25(10):1601–12.

61. Castillo JA, Epps CW, Davis AR, Cushman SA. Landscape effects on gene flow for a climate-sensitive montane species, the American pika. Molecular Ecology. 2014;23(4):843–56. pmid:24383818

62. Landguth EL, Cushman SA, Schwartz MK, McKelvey KS, Murphy M, Luikart G. Quantifying the lag time to detect barriers in landscape genetics. Molecular Ecology. 2010;19(19):4179–91. pmid:20819159

63. Shirk AJ, Cushman SA. Spatially-explicit estimation of Wright’s neighborhood size in continuous populations. Frontiers in Ecology and Evolution. 2014;2.

64. Frankham R, Briscoe DA, Ballou JD. Introduction to conservation genetics: Cambridge University Press; 2002.

65. Cheyne SM, Stark DJ, Limin SH, Macdonald DW. First estimates of population ecology and threats to Sunda clouded leopards Neofelis diardi in a peat-swamp forest, Indonesia. Endangered Species Research. 2013;22:1–9.

66. Macdonald DW, Burnham D, Hinks AE, Wrangham R. A problem shared is a problem reduced: seeking efficiency in the conservation of felids and primates. Folia primatologica; international journal of primatology. 2012;83(3–6):171–215. Epub 2013/02/01. pmid:23363584.

67. Macdonald EA, Hinks AE, Weiss D, Dickman A, Burnham D, Sandom C, et al. Identifying Ambassadors for conservation marketing. Global Ecology and Conservation. 2017;12:204–14.

68. Baguette M, Blanchet S, Legrand D, Stevens VM, Turlure C. Individual dispersal, landscape connectivity and ecological networks. Biological Reviews. 2013;88(2):310–26. pmid:23176626

69. Noon BR, McKelvey KS, Dickson BG. Multispecies conservation planning on US federal lands. Models for planning wildlife conservation in large landscapes/Millspaugh Joshua, Thompson Frank R. 2009.

70. Beier P, Majka DR, Spencer WD. Forks in the road: choices in procedures for designing wildland linkages. Conservation biology: the journal of the Society for Conservation Biology. 2008;22(4):836–51. Epub 2008/06/12. pmid:18544090.

71. Seoane J, Bustamante J, Diaz-Delgado R. Effect of Expert Opinion on the Predictive Ability of Environmental Models of Bird Distribution. Conservation Biology. 2005;19(2):512–22.

72. Shirk A, Schroeder M, Robb L, Cushman S. Empirical validation of landscape resistance models: insights from the Greater Sage-Grouse (Centrocercus urophasianus). Landscape Ecol. 2015:1–14.

73. Shirk AJ, Wallin DO, Cushman SA, Rice CG, Warheit KI. Inferring landscape effects on gene flow: a new model selection framework. Molecular Ecology. 2010;19(17):3603–19. pmid:20723066