Impact of prey occupancy and other ecological and anthropogenic factors on tiger distribution in Thailand's western forest complex

Abstract Despite conservation efforts, large mammals such as tigers (Panthera tigris) and their main prey, gaur (Bos gaurus), banteng (Bos javanicus), and sambar (Rusa unicolor), are highly threatened and declining across their entire range. The only large viable source population of tigers in mainland Southeast Asia occurs in Thailand's Western Forest Complex (WEFCOM), an approximately 19,000 km2 landscape of 17 contiguous protected areas. We used an occupancy modeling framework, which accounts for imperfect detection, to identify the factors that affect tiger distribution at the approximate scale of a female tiger's home range, 64 km2, and site use at a scale of 1‐km2. At the larger scale, we estimated the proportion of sites at WEFCOM that were occupied by tigers; at the finer scale, we identified the key variables that influence site‐use and developed a predictive distribution map. At both scales, we examined key anthropogenic and ecological factors that help explain tiger distribution and habitat use, including probabilities of gaur, banteng, and sambar occurrence from a companion study. Occupancy estimated at the 64‐km2 scale was primarily influenced by the combined presence of all three large prey species, and 37% or 5,858 km2 of the landscape was predicted to be occupied by tigers. In contrast, site use estimated at the scale of 1 km2 was most strongly influenced by the presence of sambar. By modeling occupancy while accounting for imperfect probability of detection, we established reliable benchmark data on the distribution of tigers in WEFCOM. This study also identified factors that limit tiger distributions; which managers can then target to expand tiger distribution and guide recovery elsewhere in Southeast Asia.

WEFCOM is part of an even larger conservation landscape in the Tenasserim Range, which forms Thailand's western border with Myanmar (WEFCOM, 2004). Because of its size and geographic extent, WEFCOM provides critical habitat for tigers in the region. The landscape also supports numerous other threatened wildlife species including the tiger's main prey: gaur (Bos gaurus), banteng (Bos javanicus), wild water buffalo (Bubalus arnee), and sambar (Rusa unicolor) (Simcharoen, Simcharoen, Duangchantrasiria, Bump, & Smith, 2018). Data from camera-trapping and radio telemetry suggest that HKK alone supports an estimated 35 to 58 tigers (Simcharoen et al., 2007;Duangchantrasiri et al., 2016). Female home ranges there vary from 35 to 105 km 2 . Elsewhere in WEFCOM, tiger distribution and abundance and the quality of habitat are less well known. However, ranger patrols suggest tigers are likely absent from large parts of WEFCOM, and their major prey is similarly reduced (Jornburom, 2016). This lack of knowledge restricts planning efforts to monitor and manage other protected areas in this landscape. To prioritize and strengthen the future protection of tigers throughout WEFCOM, it is important to determine where tigers occur.
Recent studies using occupancy modeling have greatly increased understanding of how habitat connectivity, human disturbance, and prey availability affect tiger occurrence. Important factors impacting tigers include forest type and extent of vegetation cover (Sunarto, Kelly, Parakkasi, & Hutajulu, 2015), connectivity of habitat (Joshi, Vaidyanathan, Mondol, Edgaonkar, & Ramakrishnan, 2013;Yumnam et al., 2014), and anthropogenic impacts such as livestock and human settlements (Harihar & Pandav, 2012;Karanth, Gopalaswamy, Kumar, Nichol, & MacKenzie, 2011;Sunarto et al., 2015). In most of WEFCOM and elsewhere in Thailand, a significant knowledge gap remains about the tiger's distribution. Previous tiger-habitat relationship studies throughout WEFCOM were limited and used inconsistent spatial and temporal scales (Kanchanasaka et al., 2010;Rabinowitz, 1993;Smith, Ahern, & McDougal, 1998;WEFCOM, 2004). They also neglected to account for imperfect detection, which can obscure the underlying ecological processes that determine distribution and habitat relationships, especially when surveys are conducted over a large landscape (Pollock et al., 2002). Therefore, it is difficult to reliably compare results across these studies. To establish reliable benchmark data on tiger distribution patterns, we conducted spatially replicated occupancy surveys throughout WEFCOM, similar to those being applied across other tiger landscapes (Harihar & Pandav, 2012;Hines et al., 2010;Karanth et al., 2011).
Occupancy models applied to multiple spatial scales, with both landscape and fine-scale predictors, have been demonstrated to be effective in addressing wildlife conservation needs (Scott et al., 2002) and have been used specifically to help inform conservation strategies of multiple species, including tigers (Wikramanayake et al., 2011). At a regional level, landscape-wide assessment of tiger distribution facilitated the identification of source populations, metapopulation structure, and functional corridors that allow individuals to move through habitat impacted by human disturbances (Karanth et al., 2011;Ranganathan, Chan, Karanth, & Smith, 2008;Smith et al., 1998). At a finer scale, space-use patterns in source areas provide insights into local factors driving habitat use, which can help inform local management options (Sunarto et al., 2015;Vinitpornsawan, 2013). In summary, analyses that include multiple spatial scales can improve understanding of tiger-habitat relationships (Johnson, 1980;McDonald, Erickson, Boyce, & Alldredge, 2012).
In this study, we seek to: (a) identify how anthropogenic pressures, landscape features, and prey occupancy determine tiger distribution and habitat use, (b) estimate the proportion of area occupied (true occupancy) by tigers at a landscape scale using a 64 km 2 grid, (c) model habitat relationships of tigers at a 1-km 2 scale to determine drivers of habitat use, and (d) develop a predictive distribution map for tigers based on spatially explicit site use models. We conclude by providing guidance on possible alternative management options.  (Figure 1). Collectively, these three large ungulates comprise 88%-95% of the tiger's prey in WEFCOM (Pakpien et al., 2017;Simcharoen et al., 2018). This landscape covers 19,600 km 2 and consists of 17 contiguous protected areas making it the largest intact protected area in Thailand and all of southern Asia (Figure 1).

| Tiger occupancy survey
Our data were collected as part of a multispecies large mammal survey that included elephants, tigers and gaur, banteng, and sambar, the main prey of tigers. We used a 64 km 2 grid, the approximate size of a female tiger home range, to determine tiger occupancy. Within each grid cell, linear survey routes, composed of 1-km spatial replicates, were delineated with the number of replicates proportional to the amount of forest habitat within the grid cell. A maximum route length of 15 km per 64-km 2 grid cell was used if the entire cell was forested and contained no villages or agricultural land. Grid cells with <10% forest cover were not surveyed under the assumption that they were unlikely to be occupied by tigers. To determine habitat use in addition to occupancy, each 1-km survey unit was divided into 100-m spatial replicates. Thus, to estimate occupancy there were a maximum of 15 1-km replicates and to estimate habitat use within 1-km 2 grid cells there were ten 100-m replicates. Detections (e.g., F I G U R E 1 Study area and designed sample units of tiger in Western Forest Complex (WEFCOM), Thailand (2010)(2011)(2012). The map shows the spatial distribution of surveyed grid cells (those with >10% forest cover). Inset: location of the study area in Thailand is outlined by a red box. NP: National Park; WS: Wildlife Sanctuary. direct sightings, scats, pugmarks/tracks, carcasses, scent marks, and vocalizations) were recorded in each 100-m subsegment on typical tiger travel routes that included wildlife trails, mineral licks, forest roads, and river banks. Additionally, in each 100-m subsegment, substrate condition, habitat type, human activities, and evidence of domestic animals were recorded. Because of the high logistic cost of traveling to random transects, we adopted the widely used modification that reduced travel time by surveying linked replicates along a linear route (Aing, Halls, Oken, Dobrow, & Fieberg, 2011;Hines et al., 2010). Surveys were conducted in the dry season (January-May, October-December) to ensure that scat persistence was consistent.
Surveying during the dry season also helps reduce heterogeneity in detection probability induced by rainfall variation (Royle & Nichols, 2003).
This sampling design allowed us to analyze occupancy at two spatial scales or second-order habitat selection and site use or thirdorder habitat selection (Johnson, 1980). For occupancy, we used 64-km 2 grid cells as "sites" and 1-km segments as spatial replicates.
For site use, 1-km 2 grid cells were used as sites and 100-m segments were used as spatial replicates.

| Selection of ecological and anthropogenic variables
We extracted Geographic Information System (GIS)-based ecological covariates from GIS public domain data and DNP's WEFCOM database. To predict occupancy and site use of tigers, we considered five groups of factors: (a) availability of prey, (b) human disturbance, (c) forest covers (d) distance to streams, and (e) terrain (for sources and further details see Supporting Information Tables S1-S13).
Prey covariates were obtained from a companion analysis that estimated probabilities of habitat use (ψ) at the site scale of (1-km 2 ), while estimating detection (p) at the scale of 0.1 km. All prey was the probability of use by at least 1 of the 3 large prey species. Further, we hypothesized that tigers would be less likely to occur in parts of WEFCOM that experience high human disturbance. We also expected human activity would have a greater impact on occupancy than correlates related to habitat heterogeneity. To assess the impact of human activities, we included a measure of relative abundance of domestic animals (domestic), distance from villages (village), and distance from roads (road). Forest areas that have been converted to shifting cultivation (agri) inside WEFCOM were considered as a measure of habitat degradation.
For habitat-related covariates, we used Thematic Mapper™ data to calculate the proportion of forest (forest), and four separate continuous variables: percentage of hill evergreen forest (HE), dry evergreen forest (DE), mixed deciduous forest (MD), and dry dipterocarp forest (DD) (WEFCOM, 2004). For geophysical variables, we used a digital elevation model to obtain elevation, slope, and terrain ruggedness. In addition, distance to rivers and streams (stream) and low slope areas (<10% slope) within 1-km or 3-km buffers along streams (flat1km or flat3km) were also included as covariates that reflect habitat quality near streams (Linkie, Chapron, Martyr, Holden, & Leader-Williams, 2006). These spatial covariates were calculated using ArcMap 10.3 (ESRI) and ERDAS IMAGINE 2013 software.
Field-based covariates were collected by surveyors as they walked transects searching for signs of humans and domestic animals including tracks, tree cutting, gun shells, and campfires.
Substrate condition (SUB), which was hypothesized to impact detection, was recorded as soft soil, hard soil, or leaf litter. The presence of domestic animals (domestic animal) was a binary variable coded as "1" for presence or "0" for absence.

| Data analyses
We applied a first-order Markovian model (Hines et al., 2010;Karanth et al., 2011) to account for spatial dependence between adjacent replicates. This model compensates for the lack of spatial independence of replicate surveys by including segment-level occupancy (parameterized by θ 0 , θ 0 and θ π ) and detection probability (p) conditional on neighboring segment-level occupancy. The spatial dependence in segment-level occupancy (θ) is captured by θ 0 if the species is present locally but was not present in the previous spatial replicate, θ 1 if the species is present locally and was present in the previous spatial replicate, and θ π for the first replicate where there is no prior information to inform segment-level occupancy.
We developed models depicting habitat relationships at two spatial scales to better understand tiger ecology and to address conservation needs (Johnson, 1980). We used different subscripts on ̂ to differentiate estimates at these two scales; ̂6 4 refers to occupancy probability at the home-range (64-km 2 ) scale, whereas ̂1 refers to probability of site use at the 1-km 2 scale. We modeled probabilities of occupancy (ψ)and detection (p) of tigers at both scales as linear functions of the above-mentioned covariates using a logit link function (MacKenzie et al., 2002).

| Model development
To analyze occupancy and site use, we used a sequence of four steps described in more detail below. First, we examined all potential covariates and eliminated highly correlated variables. Next, we grouped ecological correlates into five categories and selected one covariate in each category for subsequent analysis. These variables were additive with no interaction terms. These five variables served as a global occupancy model (Burnham & Anderson, 2002;Duren, Buler, Jones, & Williams, 2011). Third, we held the global model constant to select the most important detection covariate for each scale. Finally, we used the best detection model from step 3 to examine our candidate set of models to identify the most important predictors for occupancy and site use (Karanth et al., 2011;Royle & Nichols, 2003).
To select the appropriate covariates to develop occupancy models, we first explored the correlations among ecological variables to avoid collinearity. We considered correlation coefficients <0.7 as acceptable to include two ecological variables (Dormann et al., 2013).
We grouped variables into five categories: availability of prey, human disturbance, forest cover, proximity to stream, and terrain (Table 1). To select the factors included in our "global model," we modeled covariates in each category in a univariate fashion and chose the factor with the lowest AIC value in each category to incorporate into our global model.
We then used the global model for occupancy to evaluate variables that influenced detection. Detection probability (p) was modeled as a function of substrate type, presence of domestic animals, and area of low slope (<10% slope) within 1 km or 3 km from streams.
In this step, we held occupancy constant using the "global model" (MacKenzie et al., 2006).
The final step was to use the best model for detection to evaluate a set of models to predict occupancy at the 64 km 2 scale and site use at 1 km 2 scale (Karanth et al., 2011). For these models, we kept θ 0 , θ 1 , and θ π constant (Karanth et al., 2011;Thapa & Kelly, 2017). Estimates of coefficients and standard errors (β, SE(β)) were used to determine effect sizes and direction of influence of covariates on the probabilities of occupancy, site use, and detection. Prior to modelling, all continuous covariates were centered and scaled, (x i −x i )∕SD(x i ), to facilitate estimation of parameters using numerical optimization techniques (Donovan & Hines, 2007) and to facilitate comparisons among competing variables.
Estimates of regression parameters and the associated variancecovariance matrix derived from program PRESENCE version 12.7 (Hines, 2006) were used to explore the effect of individual covariates while holding all other covariates constant. We calculated variances of predicted values using the delta method as implemented in the CAR package in R (Fox, 2016).

| Estimation of overall tiger occupancy
To estimate overall tiger habitat occupancy (̂6 4 ) within WEFCOM (total proportion of the landscape occupied by tigers taking into account imperfect detection), we used the top ranked model (Burnham & Anderson, 2002) and averaged predictions from all 309 grid cells.
Because our sampling was a near-complete survey of the WEFCOM landscape with 309 grids of 64 km 2 , we were able to estimate overall tiger occupancy in the WEFCOM as � Ψ = ∑ 309 i=1 a îi 15,672 where a i is the forested area in cell i (i.e., total area of potential tiger habitat is 15,672 km 2 ) (Karanth et al., 2011;Srivathsa, Karanth, Jathanna, Kumar, & Karanth, 2014). We used a parametric bootstrap (Efron & Tibshirani, 1994) to compute covariance and the standard error of overall tiger occupancy (̂6 4 ).
To map tiger distribution, we employed our best-supported model for tiger site use (1 km 2 ) to help managers identify the key factors impacting spatial distribution of tigers (Lakshminarayanan, Karanth, Goswami, Vaidyanathan, & Karanth, 2015). Mapping at the scale of occupancy (64 km 2 ) or use (1 km 2 ) produced very similar maps, but the finer-scale map provides better visualization of habitat use that managers need in making decisions. We used a probability of use >0.6 as a convenient metric to indicate high-quality tiger habitat.

| RE SULTS
We surveyed a total of 3,517 1-km segments distributed in 309 (64 km 2 ) grid cells across WEFCOM to determine occupancy.

| Influence of covariates on tiger occupancy and site use
The global model at the 64-km 2 scale included the variables: all prey, proportion of forest, elevation, distance from streams, and relative abundance of domestic animals. At the 64-km 2 scale of site use, the same ecological variables were chosen, except "all prey" was replaced by "sambar." At both scales, we found that the model incorporating low slope areas along streams was the best model for detection.
At the home-range scale (64-km 2 ), all prey together (i.e., gaur, banteng, and sambar) was the most important predictor of occupancy based on the size of standardized regression coefficients (β) ( Table 1).
Other important factors for predicting occupancy were proportion of forest, elevation, relative abundance of domestic animals, and distance from streams. The AIC-best model with 58% of the model weight revealed that prey availability, proportion of forest, and elevation were positively correlated with tiger occupancy, and relative abundance of domestic animals was negatively correlated with tiger occupancy (Table 1). The probability of tiger occupancy increased from 20% to 80% as the probability of "all prey" site use increased from approximately 30% to 80% when holding other variables at their mean ( Figure 2). A complete set of 32 models and model-specific regression coefficients are presented in Supporting Information Tables S1-S13.
For site use (1 km 2 ), the most important predictor was sambar presence. In addition to sambar, the best-supported model included a negative effect of distance from streams and negative effect of domestic livestock. This model garnered 99% of the model weights (Table 1). The probability of tiger site use increased from 20% to 80% as the probability of sambar occupancy increased from approximately 40% to 75% when holding other predictors at their mean values (Figure 3). Model selection results and coefficient estimates are presented in Supporting Information Tables S1-S13.
We generated a predicted distribution map using model-based probabilities of tiger site use ≥0.60 (ranging from 0 to 1) from the best-supported model for site use. Areas of predicted tiger site use (≥0.60) were restricted to the east-central and northeastern regions of WEFCOM (Huai Kha Khaeng, east and west Thung Yai, Umpang, and Mae Wong) (Figure 4). Site use throughout much of the remainder of the landscape was less contiguous, consisting largely of scattered '"islands"' of predicted (site use) ≥0.60.

| Estimate of tiger occupancy
We estimated tiger occupancy ̂6 4 to be 0.37 (SE 0.06) for the 309 surveyed cells or 5,858 km 2 (SE 1,758 km 2 ) of the 15,600 km 2 of potential habitat. This estimate was 21% larger than the naïve estimate of 0.27.

| D ISCUSS I ON
Our study is the first modeling effort to incorporate occupancy of prey species as a covariate in tiger occupancy models. Despite widespread use of occupancy modeling, very little research has been published on how prey occupancy influences the distribution F I G U R E 2 Relationship between the highly influential covariates (based on regression coefficient (β) and 95% CI from best-supported model) and the probability of tiger occupancy in WEFCOM, Thailand (2010. Effect sizes (beta estimates) are based on standardized data while holding the other covariates at their mean values. Tick marks on the X-axis show density of data values in 64 km 2 grid cell. See Supporting Information Tables 1 and 2 for description of covariates of large carnivores (Andresen, Everatt, & Somers, 2014;Everatt, Andresen, & Somers, 2015). Incorporation of prey variables into models has often been done using indices such as overall prey biomass or density (Robinson, Bustos, & Roemer, 2014), general presence or absence of prey (Alexander, Gopalaswamy, Shi, Hughes, & Riordan, 2016;Vinitpornsawan, 2013), relative abundance indices (Chanchani, Noon, Bailey, & Warrier, 2016), or photo-trap rates (Sunarto et al., 2015). Our study considered only the probability of large prey occupancy because large prey comprises 89% of tiger prey biomass (Simcharoen et al., 2018). We also examined other natural and anthropogenic features to determine the relative influence of these correlates in shaping tiger distribution in WEFCOM.

Identification of environmental and anthropogenic factors affect-
ing the distribution of tigers not only increases our understanding of tiger occupancy, but also helps target those correlates that can be managed to increase tiger distribution (Fieberg & Johnson, 2015).
Our study analyzed occupancy at two spatial scales and yielded an important difference. At the 64-km 2 scale, the model with all large prey was the highest ranked model. Total large prey biomass was also highly inversely correlated to the size of a female tiger's home range .
However, within a tiger's home range, tigers preferred localized areas dominated by sambar. This result is consistent with the fact that across the tiger's range, sambar is both a preferred prey and the dominant prey biomass in the diet of tigers (Hayward, Jędrzejewski, & Jedrzejewska, 2012). These differences show the importance of analyzing data at multiple scales.
While tiger occupancy is largely influenced by prey availability, we also found that tiger occupancy decreased with greater distance from streams ( Figure 3). We suspect these areas are crucial for tiger site use because low slope forests near streams are also preferred habitat of sambar and banteng (Jornburom, 2016;Simcharoen, Savini, Gale, Roche, et al., 2014). Tigers also occupied areas of higher altitude F I G U R E 3 Relationship between the highly influential covariates (based on regression coefficient (β) and 95% CI from best-supported model) and the probability of tiger site use in WEFCOM, Thailand, 2010-2012. Effect sizes (beta estimates) are based on standardized data. Tick marks on the X-axis show density of data values in 1 km 2 grid cell. See Appendix 1 for description of covariates F I G U R E 4 Spatially explicit predictions map of tiger site use constructed from the best-supported occupancy model developed at 1-km 2 scale based on the analysis of occupancy surveys (2010)(2011)(2012) in WEFCOM, Thailand in the core of WEFCOM; this finding might be attributed to avoiding human activities near villages, especially in the west and the south of WEFCOM. Thus low elevation areas per se may not limit tigers, but, instead, livestock grazing at low elevation may degrade the habitat of the tiger's primary prey. Also, poaching sign is more common at lower elevations (Jornburom, 2016;Vinitpornsawan, 2013).
Our study demonstrated that large parts (63%) of the WEFCOM landscape were devoid of tigers and that tiger habitat use was concentrated in core protected areas in the north (Figure 4), whereas there were only a few scattered patches that were identified as potential habitat in southern areas where tigers were historically widely distributed (Smith et al., 1999). These current patterns of tiger distribution show a positive response to the presence of large ungulates and a negative response to domestic cattle grazing.
Previous research has shown the relationship of tiger density to prey density (Karanth, Nichols, Kumar, Link, & Hines, 2004;. Given that distribution and abundance of large ungulates are critical to the distribution and abundance of tigers, additional research is needed to identify the key ecological correlates that drive ungulate distribution and abundance in western Thailand.

| CON CLUS IONS
Our results suggest that low tiger occurrence in WEFCOM is primarily due to low abundance of large prey and Jornburom (2016) attributed their low abundance to human activities near villages.
Several other occupancy studies note that scarcity of natural prey near villages is a consequence of degradation of habitat by livestock grazing (Harihar & Pandav, 2012;Karanth et al., 2011). However, relocating villagers, who have long historical residence, may not be an accepted management strategy. Furthermore, globally there can be strong opposition to forced resettlement (Clements, Suon, Wilkie, & Milner-Gulland, 2014. The well-established Smart Patrolling system has had an overall strong protection impact in WEFCOM, but has not increased prey abundance or eliminated domestic livestock grazing or subsistence poaching near villages (Duangchantrasiri et al., 2016). Thus increasing tiger population size may depend on reducing certain activities in the vicinity of villages.
Jornburom (2016)  Thung Yai West (two villages) with the goal of establishing co-management (United Nations Development Program, 2015). Specific objectives of co-management in these and future studies could be beneficial to tiger conservation. In a similar situation in Nepal, co-management has been effective in reducing poaching and habitat degradation in forests in close proximity to villages (Carter, Shrestha, Karki, Pradhan, & Liu, 2012).
Our research shows that the low occupancy of tigers (37%) was largely a consequence of the absence of large prey. Managers need to identify the factors that limit the distribution of large mammals and test new options to increase large ungulate distribution so that WEFCOM remains source population able to support tiger restoration efforts elsewhere in Thailand and Southeast Asia.

ACK N OWLED G M ENTS
We

DATA ACCE SS I B I LIT Y
Presence/Absence of tiger in WEFCOM Thailand: Dryad https://doi.