Assessing the adequacy of a protected area network in conserving a wide‐ranging apex predator: The case for tiger (Panthera tigris) conservation in Bhutan

Protected area networks (PAN) are essential for conserving wide‐ranging apex predators but their adequacy in species protection has rarely been assessed. Here, we assess the adequacy of Bhutan's PAN in conserving and providing connectivity to the endangered tiger (Panthera tigris). We determine the current extent of tiger habitat, predict new suitable habitat, identify potential corridors, and empirically estimate the range of tiger numbers that Bhutan can spatially support. We use two spatial models with different approaches to ascertain current tiger distribution and predict new suitable tiger areas: (a) an expert model based on tiger ecology and (b) an observation model from observed tiger distribution. The expert model identified more suitable tiger areas (32,887 km2) over the observation model (29,962 km2), with the PAN encompassing 46% and 45% of predicted suitable areas, respectively. Vast suitable tiger habitat remains unprotected. Based on our estimates of total suitable habitats, Bhutan can spatially support 138–151 tigers compared to the current estimate of 103, thereby precluding a doubling in tiger numbers. To ensure adequate protection of tigers in Bhutan, we recommend readjusting and/or expanding existing PAN boundaries, including the designation of new corridors, protecting habitats, and conserving prey populations.

Species distribution and or habitat suitability models have identified suitable habitats for priority species (Raxworthy, Martinez-Meyer, Horning, et al., 2003), and they are useful tools for species conservation planning (Chefaoui, Hortal, & Lobo, 2005;Hirzel, Le Lay, Helfer, Randin, & Guisan, 2006). The last two decades has seen an increase in GIS-based predictive models (species distribution models, habitat suitability models, and niche models) which assign suitability values (Hirzel & Le Lay, 2008) to geographic areas for spatially predicting a species' presence/absence. These models subsequently identify suitable areas for conserving priority species amidst rapidly changing landscapes . Although ascertained as reliable (Hirzel et al., 2006), questions remain on the veracity of identified suitable area and whether they adequately reflect the full breadth of a species' niche to meet fundamental ecological requirements.
Here, we highlight the case of the endangered tiger (Panthera tigris), Asia's apex predator and flagship species (Wang, Royle, Smith, et al., 2018). It has experienced substantial population decline and range contraction from habitat destruction and fragmentation, illegal trade in body parts, and conflict with humans (Wikramanayake, Dinerstein, Seidensticker, et al., 2011), despite providing ecological benefit through reducing livestock and crop depredation (Thinley, Rajaratnam, Lassoie, et al., 2018). Tiger populations have dramatically reduced from close to 100,000 in the early 20th century to less than 3,600 in the early 21st Century, prompting a global conservation initiative to double tiger numbers by 2022 (GTIS, 2011). Their population viability is at a critical point and conservation strategies implemented over the next decade will determine their fate beyond existence in protected areas (Chanchani et al., 2016;Dutta et al., 2016), because their metapopulation dynamics are subject to effective dispersal of individuals between habitat patches in the landscape mosaic Wang et al., 2018). Therefore, it is crucial to not only maintain tiger populations at key sites, but also enable their persistence across much larger landscapes through well-linked habitat networks where survival and reproduction are complemented by opportunities for dispersal and colonization (Gubbi, Mukherjee, Swaminath, & Poornesha, 2016;Thinley et al., 2020).
Because of its status as a national priority species (Wang & Macdonald, 2009), landscape-level tiger studies in Bhutan have mainly focused on population (Tempa, Hebblewhite, Goldberg, et al., 2019;Thinley, Dorji, Tempa, et al., 2015) and movement dynamics , as conducted elsewhere in South India (Gubbi et al., 2016), China (Wang et al., 2018), and the Terai Arc Landscape straddling Nepal and India (Thapa, Wikramanayake, Malla, et al., 2017). In an effort to maintain viable tiger populations, some studies have additionally prioritized assessment of landscape connectivity through various approaches such as mapping potential forest loss in Sumatra (Poor, Shao, & Kelly, 2019), investigating metapopulation gene flow in central India (Seidensticker, 2016), and spatially identifying corridors in western India (Mondal, Habib, Talukdar, & Nigam, 2016). More recently, Sanderson, Moy, Rose, et al. (2019) even investigated shared socioeconomic pathways to gauge the influence of human population growth on tiger conservation.
Various modeling approaches have been increasingly utilized to promote effective tiger conservation. These include assessing connectivity for tigers through leastcost corridor modeling and circuit theory in central India (Dutta et al., 2016), investigating limiting effects of landscape features on tiger gene flow in the Western Ghats of India (Reddy, Puyravaud, Cushman, & Segu, 2019), determining population viability through genetic modeling in central India (Thatte, Joshi, Vaidyanathan, Landguth, & Ramakrishnan, 2018), and occupancy modeling based on anthropogenic and environmental covariates to estimate occupancy and habitat use by tigers in China (Wang et al., 2018) and the Central Terai Landscape of India (Chanchani et al., 2016). Modeling in Bhutan for tiger conservation has been limited. Rostro-García et al. (2016) modeled predation risk to mitigate livestock predation by tigers (Rajaratnam, Vernes, & Sangay, 2016). Recently, Tempa et al. (2019) used a Bayesian-based spatially explicit capture-recapture model to produce a tiger density map, while Penjor, Tan, Wangdi, and Macdonald (2019) used hierarchical occupancy modeling to understand the environmental and anthropogenic correlates of tiger presence. Despite Bhutan harboring a healthy population of tigers within a well-connected PAN (Thinley et al., 2015), there remains the question of protected area effectiveness (Hausner, Engen, Bludd, & Yoccoz, 2017) to meet conservation targets specific to tiger conservation.
The principal aim of our study is to assess the adequacy of Bhutan's PAN for tiger protection and connectivity. To achieve this, we model the extent of tiger habitat and assess their coverage in Bhutan's PAN. Additionally, we compare the extent of tiger habitats yielded by our predictive models to ascertain new suitable tiger areas, identify new potential corridors for effective tiger conservation, and empirically estimate the range of tiger numbers Bhutan can spatially support.

| Study area
The Eastern Himalayan Kingdom of Bhutan ( Figure 1) is synonymous with strong conservation policy and biodiversity richness. Within a small geographic size of 38,394 km 2 , Bhutan maintains a contiguous forest cover of approximately 71% (Thinley, Rajaratnam, Tighe, et al., 2019), which supports a rich biodiversity of 11,248 species (NBC, 2017). People, additionally, revere nature from the influence of Buddhist philosophy (Thinley et al., 2019). Because Bhutan's Constitution mandates a permanent forest cover of at least 60% of total land area, approximately 51% of the country is protected : 41% as protected areas and 10% as biological corridors for wide-ranging species, such as tigers, common leopards (Panthera pardus), elephants (Elephas maximus), and dholes (Cuon alpinus; Figure 1). With an estimated tiger population of 103 individuals at a density of 0.46 individuals per 100 km 2 (Thinley et al., 2015), Bhutan is an acknowledged stronghold for tigers in the Eastern Himalayas (Tempa et al., 2019).

| Modeling suitable areas for tiger
We developed two predictive models on tiger distribution: a rule-based model (hereafter called the expert model) based on expert knowledge of tiger ecology (Johnsingh, 2004;Karanth, 2001;Schaller, 1967;Seidensticker, 1996) without utilizing tiger occurrence points or its current distribution; and an observationbased model (hereafter called the observation model) using correlations between observed tiger occurrence and environmental variables. The expert model was based on an Analytic Hierarchy Process (AHP) framework (Saaty, 1987) which assigned weights derived by a square matrix with values inserted from expert judgment, to spatial layers in ArcGIS. AHP-based modeling is widely used in multi-criteria decision making, resource planning, and fire hazard mapping (Chhetri & Kayastha, 2015 Thinley, Shafapour, Thinley, & Shabani, 2020). It is a decision-making tool which compares and weights multiple criteria or factors influencing a certain phenomenon. In modeling our tiger distribution, the AHP matrix provided a logical step in determining influential weights of each spatial layer which is comparable to conventional rule-based modeling (e.g., Schadt et al., 2002) where experts assign weights to the layers. The observation model was based on presence-only records (Phillips, Anderson, & Schapire, 2006) of tiger occurrence, and constructed using MaxEnt program (Version 3.4.1) and ArcGIS. MaxEnt is a species distribution modeling program that uses presence-only data (georeferenced presence points) of a species along with environmental covariates (selected spatial layers in our case) to generate a spatial surface of probability of occurrence or suitability (Elith et al., 2010;Phillips et al., 2006). Based on tiger ecology, anthropogenic influence on tigers, and availability of spatial layers, we assigned eight key input variables to both models. Variables were processed in ArcMap version 10.7.1 using a standard cell size of 30 m × 30 m corresponding to the lowest resolution of Bhutan's Digital Elevation Model (DEM) as follows: a. "Elevation": prepared by reclassifying Bhutan's DEM (Jarvis, Reuter, Nelson, & Guevara, 2006) into ten categories ranked for tiger suitability based on accessibility, climatic regime, and temperature. Lower elevations were ranked more suitable than higher elevations (Table S1). b. "Slope": extracted from the DEM and classified into ten categories based on their physical influence on tiger movement. Flat terrain was ranked more suitable than steep terrain due to energetic cost associated with negotiating steep terrain (Table S1), as evidenced by negative correlation between tiger gene flow and rugged terrain in India's Western Ghats (Reddy et al., 2019). c. "Prey": prepared by merging the distribution of nine prey species of sambar (Rusa unicolor), gaur (Bos gaurus), muntjac (Muntaicus muntjac), water buffalo (Bubalus arnee), goral (Naemorhedus goral), serow (Capricornis thar), hog deer (Axis porcinus), musk deer (Moschus chrysogaster), and takin (Budorcas taxicolor), and reclassified as a single layer comprising 10 categories. We used individual prey distribution maps from the Field Guide to the Mammals of Bhutan (Wangchuk, Thinley, Tshering, et al., 2004) which were derived from presence records obtained through field surveys and museum collections. Higher values were assigned to areas with higher prey diversity (Table S1), because prey abundance is an important ecological determinant of tiger density in tiger landscapes (Karanth, Nichols, Kumar, Link, & Hines, 2004;Karanth & Sunquist, 1995). d. "Drainage," that is, distance from rivers and streams: sourced from the Drainage Map of Bhutan and rasterized using the Euclidean Distance tool. Higher values were assigned to areas closer to rivers and streams which are frequented by tigers (Table S1). e. "Settlement," that is, distance from human settlements: created from the Settlement Map of Bhutan 2006 (OCC, 2005) and rasterized using the Euclidean Distance tool. Higher values were assigned for areas situated farther away from settlements (Table S2) because tiger presence is positively associated with greater distance from human settlements (Penjor et al., 2019). f. "Land-use": derived from the Land-use Map of Bhutan 2011 (NSSC & PPD, 2011) encompassing major forest types, agricultural lands, and other land-use types.
Higher ranking was assigned to land-use features deemed more suitable to tigers (Table S2), for example, densely forested areas which provide cover to ambush prey (Seidensticker, 1996). Conversely, lower ranking was assigned to relatively open pine forests which are rarely frequented by tigers (Schaller, 1967). g. "Road," that is, distance from roads: prepared by rasterizing the latest Road Map of Bhutan obtained from the Ministry of Works and Human Settlement in a similar approach used for the settlement layer, with higher suitability scores for areas occurring farther away from roads (Table S2) due to negative correlation between tiger occurrence and distance to roads (Linkie, Chapron, Martyr, Holden, & Leader-Williams, 2006). h. "Religion": Euclidian distance from religious sites (Buddhist temples, monasteries, citadels of deities, and other religiously significant areas) digitized from Google Earth™. Higher values were assigned to areas situated closer to religious sites due to strong Buddhist sentiment toward all sentient beings (Table S2), as evidenced by snow leopard affinity to Buddhist monasteries in Tibet (Li, Wang, Yin, et al., 2013).
For the expert model, we assigned paired comparison weights to variables in an 8 × 8 AHP matrix (Table 1) based on expert judgment (input from local tiger experts and experienced forestry professionals) and knowledge of tiger ecology. Weights were recalibrated until a consistency ratio (CR) of lesser than 0.10 (0.04 in our case; Table 1) was obtained after several iterative processes in divergence from a principal eigenvalue. These weights were then entered as percentages of influence from each spatial layer in ArcGIS (Tshering et al. 2020) to generate a tiger suitability map. The AHP matrix assigned the highest weight to the prey layer followed by the Landcover layer (Table 1), which was consistent with tiger ecology since prey abundance is the primary determinant of tiger survival (Karanth & Stith, 1999) while cover is important for ambushing prey (Schaller, 1967;Seidensticker, 1996). For the observation model, we utilized a database of 145 georeferenced tiger presence points in Bhutan from 2014 to 2019 (Figure 1). Points were derived from a 2014 to 2015 comprehensive nationwide camera trapping tiger survey (Thinley et al., 2015), ancillary park-wide surveys in intervening years, and adhoc field reports of tiger evidence (indirect signs, livestock kills, and direct sightings). We omitted two points occuring with others within the same cell size of 30 meters to minimize errors resulting from spatial autocorrelation (Kanagaraj, Wiegand, Kramer-Schadt, Anwar, & Goyal, 2011). We then used the tiger occurrence points in MaxEnt by correlating them with input variables (Table 2) to produce a probability surface of tiger occurrence (Elith et al., 2010). In line with Namgyal and Thinley (2017), we used the default settings in the MaxEnt model comprising 500 interactions, a 0.00001 convergence threshold, 1 regularization multiplier, and 10,000 background points with 50% random test percentage. This was exported to ArcMap and reclassified into a suitability surface. Both models generated maps with continuous suitability gradients (Hirzel et al., 2004). They were normalized (with probability values ranging from 0 to 1) and classified into three classes ("highly suitable," "moderately suitable," and "unsuitable") using the natural breaks in ArcMap. Consequently, values from 0 to 0.1 were classified as unsuitable; 0.1 to 0.5 as moderately suitable; and >0.5 as highly suitable.
Using a three-way error/confusion matrix, we assessed the expert model through actual field verification of suitability classes within a 30 m × 30 m area at 500 random locations in each suitability class. We rapidly assessed tiger suitability in the field using the following criteria: overlap/presence of human structures, rocky outcrops, mining areas, slope >75 degree, and absence of prey evidence (direct and indirect signs) as "unsuitable"; absence of man-made structures, presence of pure conifer forests, plantations, pastures, and agriculture, slope between 35 and 75 degrees, and presence of only one prey species as "moderately" suitable; and absence of human structures, presence of broadleaved, mixed conifer, and sub-alpine forests, slope < 35 degrees, and presence of ≥ 2 prey species as "highly suitable." The observation model was assessed using the area under the receiver-operating curve (AUC; Pourghasemi, Pradhan, & T A B L E 1 An 8 × 8 Analytic Hierarchy Process matrix showing pairwise comparison of eight variables (LC, Land cover; RS, distance from religious sites; RD, distance from road; HS, distance from human settlement; WB, distance from water bodies; PD, prey diversity; EL, elevation; SL, slope) with respect to perceived influence on tiger (Panthera tigris) niche in Bhutan, based on expert knowledge Gokceoglu, 2012) using 50% of the tiger observation points as training data and the remaining 50% as test data for validation. AUC values from 0 to 0.5 indicate poor fit of the data while values closer to 1 indicate perfect fit (Fielding & Bell, 1997).

| Assessing adequacy of tiger protection and connectivity
We first computed areas (km 2 ) for each suitability class in both models. We then overlaid a layer of protected areas and biological corridors over suitability layers from both models. Adequacy of tiger protection was assessed by calculating proportion of highly suitable (optimal) and moderately suitable (suboptimal) areas captured inside protected areas. Likewise, adequacy of tiger connectivity was assessed based on the proportion of suitable areas (optimal and suboptimal combined) captured in corridors.

| Adequacy of tiger protection and connectivity
Our expert and observation models revealed that protected areas captured only 35% (11,488 km 2 ) and 33% (9,739 km 2 ) of total suitable areas for tigers, respectively (Table 3). Within these areas, expert and observation models predicted 31% (5,272 km 2 ) and 32% (3,991 km 2 ) as optimal areas for tigers, respectively. Expert and observation models additionally assessed biological corridors to encompass only 11% (3,717 km 2 ) and 12% (3,688 km 2 ) of suitable areas for tigers, respectively (Table 3). Within these areas, expert and observation models further identified 12% (2,065 km 2 ) and 15% (1,894 km 2 ) as optimal tiger areas, respectively. Contrastingly, suitable tiger areas equating to 46% (17,682 km 2 ) from the expert model and 43% (16,535 km 2 ) from the observation model, existed outside protected areas and biological corridors (Table 3). Our expert and observation models additionally deemed 53% (9,425 km 2 ) and 40% (6,560 km 2 ) of these areas as optimal for tigers.

| DISCUSSION
We highlight differences in modeling approach and outputs between our spatial predictive models. The expert model, with no reliance on current species distribution to map potential tiger niche but based on expert knowledge of tiger ecology, identified more suitable and optimal areas for tigers over the observation model for which current species distribution was integral. This highlights a potential shortcoming in ecological modeling and species conservation approaches specifically relying on current species distribution. It is possible that many tiger individuals or populations are narrowly persisting in suboptimal areas which are protected based on the current observed presence. Thus, we posit that current tiger distribution may not reflect its best needs, but rather, reflects remnants of the historic range to which the species has been relegated due to anthropogenic activities (Kanagaraj et al., 2011;Seidensticker, 2010). Additionally, observation models may often designate suitable areas based on detection of transient movement (Hirzel et al., 2004) by a wide-ranging species such as the tiger, when in fact, these areas may be deemed unsuitable by the expert model due to not meeting essential niche requirements (Hirzel & Le Lay, 2008). As such, models using only current or observed distribution can reveal environments where species can exist, but not necessarily where they can optimally thrive. Expert models guided by ecological knowledge on target species may allow sufficient assessment on whether current species distribution reflects the full breath of its niche. We, therefore, encourage ecologists and conservationists to use expert models or historical records as additional frames of reference to compare with observation models, when inferring total estimated suitable area for a species. Identification of new suitable and optimal areas overlaid with PAN maps allows us to pinpoint some conservation mismatch (Kukkonen & Tammi, 2019). In Bhutan's case, our predictive models revealed vast unprotected suitable and optimal tiger habitats. This is not surprising given the 2014-2015 National Tiger Survey of Bhutan identified more tigers outside the PAN (Thinley et al., 2015) while Tempa et al. (2019) identified some tiger source (or high density) sites outside the PAN. Similarly, Kanagaraj et al. (2011) predicted 24% (18,500 km 2 ) of the Terrai Arc Landscape of India and Nepal (a potential Himalayan tiger refuge adjoining lower Nepal and northern India) as suitable for tiger but only 7% was protected. A similar situation could hold true in other tiger range countries. As such, we demonstrate that protection and connectivity are often inadequate to support a functional metapopulation for sufficient species dispersal (Bargelt, Fortin, & Murray, 2020).
Our model outputs promote proactive tiger conservation in Bhutan. We recommend expanding or readjusting boundaries of protected areas and enhancing connectivity to encompass the full breadth of the tiger's niche. Such measures should focus on suitable and optimal tiger habitats within warm broadleaved forests in the southern foothills. Specifically, boundaries of Phibsoo Wildlife Sanctuary (PWS) and Jhomotshangkha Wildlife Sanctuary (JWS; Figure 1) need expansion to encompass more tiger habitats. A new corridor is essential to connect Bumdeling Wildlife Sanctuary (BWS) with Sakteng Wildlife Sanctuary (SWS; Figure 1), as further evidenced by documented tiger reappearance in BWS after a likely absence of 12 years . SWS is also needs to be linked to JWS with both sanctuaries recording tigers during the 2014-2015 National Tiger Survey (Thinley et al., 2015), thereby potentially harboring important metapopulations. Corridors should also be designed to connect PWS with Jigme Khesar Strict Nature Reserve, and Royal Manas National Park with SWS ( Figure 1). Current corridors connecting PWS and Jigme Dorji National Park with Jigme Singye Wangchuck National Park (Figure 1) also needs readjustment for better connectivity between tiger metapopulations and suitable habitats.
We further highlight the potential to expand tiger conservation efforts beyond protected areas (Jonas, Barbuto, Jonas, Kothari, & Nelson, 2014) in tiger range countries. In Bhutan's case, tiger habitats outside current protected areas are managed as state reserve forests by Territorial Forest Divisions (TFD) which are administrative sub-branches under the Department of Forests and Park Services which manage protected areas. TFDs which also manage biological corridors can thus be assigned new roles to protect tigers and their habitats with commensurate funding support, thereby increasing tiger protection.
Tigers lost four subspecies and 93% of their historical range (Thatte et al., 2018), and the global population plummeted to less than 3,600 individuals a decade ago (GTIS, 2011) from habitat loss, intense poaching, and inadequate conservation effort (Dinerstein, Loucks, Wikramanayake, et al., 2007). This prompted proposed conservation measures such as the "six percent solution" to protect 42 key tiger source sites encompassing 6% of tiger distribution (Walston, Robinson, Bennett, et al., 2010), and globally doubling tiger numbers by 2022 through increased conservation efforts in prioritized tiger landscapes (GTIS, 2011;Wikramanayake et al., 2011). Despite Bhutan being acknowledged as a global tiger hotspot (Tempa et al., 2019), doubling tiger numbers in Bhutan is not ecologically feasible because total suitable habitat predicted by our models can only support 138 to 151 individuals relative to the 103 tigers currently identified with certainty (Thinley et al., 2015). It is, however, a small step toward the global call to double tiger numbers by 2022. In an attempt to assess the reality of doubling tiger numbers, Harihar, Chanchani, Borah, et al. (2018) recently reviewed the recovery potential of 18 tiger source sites identified under the World Wide Fund for Nature's (WWF) Tigers Alive Initiative. They estimated that these sites currently support 165 (118-277) tigers with the potential to harbor 585 (454-739) individuals, constituting a modest 15% increase in the global population. Given the lack of enough tiger habitat to reach the goal of doubling tiger numbers in Bhutan, additional habitat protection and restoration is needed to conserve prey populations and potentially expand tiger distribution to ensure the long-term survival of this iconic species in Bhutan.