On the tiger trails: Leopard occupancy decline and leopard interaction with tigers in the forested habitat across the Terai Arc Landscape of Nepal

Better conservation planning requires updated information about leopard distribution to prioritize and allocate limited resources available. The long-term persistence of leopards and sympatric tigers can be compromised by linear infrastructure development such as roads that fragment habitat. We used detection and non-detection data collected along walking search paths (~4140 km) in 96 grid cells (each cell 15 km by 15 km) spread across potential habitat (~13,845 km2) in the Terai Arc Landscape, Nepal. Multi-season occupancy models allowed us to make both spatial and temporal inferences between two surveys in 2009 and 2013, based on ecologically relevant covariates recorded in the field or remotely sensed. Additionally, we used 2013 data to make inferences on co-occurrence between tigers and leopards at the landscape level. We found the additive model containing deforestation and district roads negatively influenced leopard detection across the landscape. Although weak, we found anthropogenic factors such as extent of deforestation (decrease in forest cover) negatively affected leopard occupancy. Road abundance, especially for the east-west highway and district roads, also negatively (but weakly) influenced leopard occupancy. We found substantially lower occupancy in the year 2013 (0.59 (SE 0.06)) than in 2009 (0.86 (SE 0.04)). Tigers and leopards co-occurred across the landscape based on the species interaction factor (SIF) estimated at 1.47 (0.13) but the amount of available habitat and the prey index mediated co-occurrence. The SIF decreased as habitat availability increased, reaching independence at large habitat patches, but leopard occupancy declined in sites with tigers, primarily in large patches. The prey index was substantially lower outside of protected areas and leopards and tigers co-occurred more strongly in small patches and at low prey indices, indicating potential attraction to the same areas when prey is scarce. Mitigation measures should focus on preventing loss of


a b s t r a c t
Better conservation planning requires updated information about leopard distribution to prioritize and allocate limited resources available. The long-term persistence of leopards and sympatric tigers can be compromised by linear infrastructure development such as roads that fragment habitat. We used detection and non-detection data collected along walking search paths (~4140 km) in 96 grid cells (each cell 15 km by 15 km) spread across potential habitat (~13,845 km 2 ) in the Terai Arc Landscape, Nepal. Multi-season occupancy models allowed us to make both spatial and temporal inferences between two surveys in 2009 and 2013, based on ecologically relevant covariates recorded in the field or remotely sensed. Additionally, we used 2013 data to make inferences on co-occurrence between tigers and leopards at the landscape level. We found the additive model containing deforestation and district roads negatively influenced leopard detection across the landscape. Although weak, we found anthropogenic factors such as extent of deforestation (decrease in forest cover) negatively affected leopard occupancy. Road abundance, especially for the east-west highway and district roads, also negatively (but weakly) influenced leopard occupancy. We found substantially lower occupancy in the year 2013 (0.59 (SE 0.06)) than in 2009 (0.86 (SE 0.04)). Tigers and leopards co-occurred across the landscape based on the species interaction factor (SIF) estimated at 1.47 (0.13) but the amount of available habitat and the prey index mediated co-occurrence. The SIF decreased as habitat availability increased, reaching independence at large habitat patches, but leopard occupancy declined in sites with tigers, primarily in large patches. The prey index was substantially lower outside of protected areas and leopards and tigers co-occurred more strongly in small patches and at low prey indices, indicating potential attraction to the same areas when prey is scarce. Mitigation measures should focus on preventing loss of

Introduction
There is increasing evidence that large carnivore populations can thrive in multi-use landscapes (Banerjee et al., 2013;Odden et al., 2014). The common leopard (Panthera pardus) is a widely distributed felid (Sunquist and Sunquist, 2002) due to its diet flexibility (Hayward et al., 2006) and its ability to live in environments ranging from rainforests to deserts, and even in areas of high human impact (Athreya et al., 2013;Gubbi et al., 2020). However, the leopard's range has collapsed considerably, now occupying only 25e37% of its historic range (Jacobson et al., 2016). The primary drivers of decline are habitat loss and fragmentation, depletion of natural prey, and direct persecution due to conflict with people, poaching, and indiscriminate killings (Athreya et al., 2011;Swanepoel et al., 2014;Rostro-García et al., 2016;Paudel et al., 2020). According to the Convention on International Trade in Endangered Species (CITES), the leopard has been identified as one of the major traded carnivores of the Asian big cats and recently, global leopard status was re-assessed and changed from Near Threatened to Vulnerable in 2016 .
The leopard (Panthera pardus fusca) and the tiger (Panthera tigris tigris) are the dominant predators in the lowlands of Nepal . They share habitat and influence trophic community structure in terrestrial ecosystems in Nepal ; India (Harihar et al., 2011); Bhutan (Wang and Macdonald, 2009); and Thailand (Simcharoen et al., 2018). The Terai Arc Landscape (hereafter Terai landscape or TAL) forms a terrestrial network of protected areas with varying degrees of forest and river connectivity among them . Multiple studies suggest highly variable leopard densities across these protected areas Thapa et al., 2014;Thapa and Kelly, 2017), but there is limited knowledge of leopard spatial distribution across the 13,845 km 2 mosaic of forests and grasslands. Information on leopard distribution and its determinants is critical for setting conservation priorities for the species (Scott et al., 2002;Jathanna et al., 2015). The leopard typifies this critical requirement for site-specific conservation planning, to adequately prioritize and allocate resources for human-wildlife conflict resolution and habitat protection measures across human-dominated areas.
Linear infrastructure development, such as the construction of roads and canals, is an emerging threat that triggers largescale habitat fragmentation and is proliferating rapidly worldwide (Trombulak and Frissell, 2000). These linear infrastructures, however, are ubiquitous features of modern landscapes and are a significant threat as they can hinder dispersal of wildlife between habitat patches surrounding urban landscapes (Forman and Alexander, 1998). Road networks can fragment large areas of continuous habitat into smaller patches (Saunders et al., 2002;Mutter et al., 2015). For carnivores in general, population densities are often low while home range sizes and dispersal distances are large (Duangchantrasiri et al., 2016;Rostro-García et al., 2016), thus carnivore persistence can be strongly influenced by the isolating effects of such linear infrastructures. Roads also can have a major effect on large-carnivore mortality directly via vehicle collisions, and through increased hunting and poaching access, and indirectly by providing greater prey hunting access to people that can result in reduced prey availability for carnivores (Kerley et al., 2002). Recent studies have reported the potential impacts of large infrastructures, particularly road development, on persistence of large mammals across the landscapes such as tigers (Carter et al., 2020) and elephants (Wadey et al., 2018). However, few studies quantified the potential impact of roads on ecological parameters such as photographic capture rates (Gubbi et al., 2012), occupancy (Chen and Koprowski, 2015), density and abundance (Fahrig and Rytwinski, 2009), and survival (Basille et al., 2013). The Government of Nepal is planning to invest US $4.7 billion in large infrastructure projects by 2030 (GoN, 2019). The impacts of these projects on carnivores need to be understood and mitigated appropriately. Regular monitoring of species is key to make timely decisions relevant to policymakers for leopard conservation in Nepal. Here, we evaluate the impacts of an existing road network on occupancy of leopards across the TAL, Nepal.
We use an occupancy modeling framework (MacKenzie et al., 2002) for quantifying distribution of leopards across the landscape. Occupancy analyses are evolving rapidly as species distribution models that consider both detection and nondetection data to explicitly account for imperfect detectability while incorporating habitat or landscape predictor variables and including other observation processes (such as surveying in potentially spatially auto-correlated trails). We employed the Hines et al. (2014) multi-season occupancy model, including spatial autocorrelation, for drawing inferences regarding occupancy dynamics of the leopard across 700 km of the flat plains mosaic of the TAL, Nepal. This modeling approach relies on both spatial and temporal replication for understanding leopard distribution over time across space. We employed various environmental, ecological, and anthropogenic variables, including those related to roads, to determine the correlates of leopard habitat occupancy along the trails originally designed for tiger surveys across the TAL, Nepal.
In addition to ecological and anthropogenic factors influencing leopard occupancy, larger predators can adversely affect guild members (Harihar et al., 2011), impacting their conservation status (Hayward and Kerley, 2008). For example, tigers are socially dominant over all co-predators including leopards, with competitive interactions affecting the time of activity, habitat use, and abundance of other sympatric carnivores (Harihar et al., 2011;Kumar et al., 2019;Lamichhane et al., 2019). Often interactions between sympatric species have been assessed at a fine-scale of a single study site, but here, we used two species occupancy models (MacKenzie et al., 2004) that included important variables from multi-season models, to assess the influence of habitat and landscape features, and species interactions between the tiger and leopard across the entire TAL, Nepal.

Study area
We conducted our study on the Nepal side of the transboundary TAL 24,710 km 2 ) stretching between the Bagmati River in the east to Shuklaphanta National Park (SuNP) in the west (Fig. 1). The landscape spans 700 linear km east to west linking protected areas along the Nepal-India border through biological corridors and the Churia hill forests ( Fig. 1)  . The TAL forms part of the Terai Duar Savannah and Grasslands ecoregion with subtropical deciduous vegetation ranging from early successional floodplain communities to climax Sal (Shorea robusta) forests (Wikramanayake et al., 2010). The TAL, Nepal contains the last remaining alluvial floodplain grasslands, which holds the highest density of tigers (0.65e16 tigers per 100 km 2 ), and leopards (3.3e5.1 leopards per 100 km 2 (Thapa and Kelly, 2017)). Besides tigers and leopards, other top predators found in the landscape include striped hyenas (Hyaena hyaena) and wild dogs (Cuon alpinus). Sambar (Rusa unicolor), swamp deer (Rucervus duvaucelii), chital (Axis axis), barking deer (Muntiacus muntjak), hog deer (Axis porcinus), wild boar (Sus scrofa), along with primates, grey langurs (Semnopithecus entellus) and rhesus monkeys (Macca mulatta), form the main prey base for the predators. Mega herbivores such as greater one-horned rhinoceroses (Rhinoceros unicornis) and Asian elephants (Elephas maximus) share the habitat with leopards in the Terai landscape. Domestic livestock, primarily cattle and goats, is found grazing in the nearby community forests and open grazing areas, and also can be potential prey for carnivores (Madhusudan, 2003). A total of 6.7 million people of mixed ethnic groups reside within 8722 km 2 of agriculture and settlement areas along the TAL, Nepal, with an annual growth rate of 1.72% (CBS, 2011). Immigration, particularly from adjoining hill districts continues to be a leading cause of human population growth in the TAL. Expanding human populations, increased development, encroachment, and extraction of firewood and fodder are the principal causes of disturbance and degradation within forests (Wikramanayake et al., 2010).

Field sampling
The occupancy survey was part of the government of Nepal's flagship monitoring protocol designed to estimate occupancy of tigers and other large carnivores at the landscape level every five years (DNPWC, 2008). Thus, while the study was originally designed for tigers (Barber-Meyer et al., 2013), we collected data along these tiger trails on multiple carnivore species across the landscape and were able to apply occupancy modeling to leopards due to the prevalence of their sign.
We divided the entire study area into 15 by 15 km grid cells (n ¼ 108). We aimed to estimate true leopard occupancy at the landscape level, thus our chosen grid cell size (225 km 2 ) was larger than a male tiger's home-range size, and also was much larger than the home range size of leopards (6e90 km 2 ) in similar habitat (Seidensticker, 1976;Norton, 1987;Simcharoen et al., 2008;Odden and Wegge, 2009). Out of the 108 grid cells, 12 were within the mid-hills and high mountain range, and these were omitted due to difficulty in surveying the undulating and steep rough terrain. Of the 96 surveyed grid cells, 32 were in protected areas, while the remaining (n ¼ 64) lay outside protected areas. The 96 surveyed grid cells covered 13,845 km 2 and included forest and grassland habitat types across the TAL, Nepal (Fig. 1). Surveys were carried out in the winter season between January and March in 2009 and 2013. Within each 15 Â 15 km grid cell, identified as site, we surveyed a maximum of 40 km of transects (defined here as search paths) and each 1 km search path was considered a spatial replicate. Thus, encounter occasions were represented by a maximum of 40 "1-km" spatial replicates. The number of 1-km replicates per grid cell was proportional to the percentage of the tiger habitat available within each cell (Barber-Meyer et al., 2013). Hence, if a grid cell consisted of only 10% tiger habitat, we surveyed at least four, 1 km search path in that grid cell. However, occasionally the realized survey effort differed from expected survey effort in the field due to logistical constraints (Harihar and Pandav, 2012).
We walked search paths within each grid cell and coded a "1" for detection and a "0" for non-detection of leopard sign. We surveyed the most strategic locations (fire-lines, roads, trails, ridges, stream beds, and riverbanks) within each grid cell to increase chances of detecting leopard and tiger signs (scats, tracks, calls, direct sightings, and kills). While there is the possibility of confusion in IDs of some leopard and/or tiger signs, we think this is low due to the much larger size of tigers (65e306 kgs) compared to leopards (28e90 kgs) (Seidensticker, 1976). Leopard signs were discerned from tiger signs based on their morphological characteristics such as much smaller sizes of pugmarks: front foot width (~average 8e9 cm for leopard, 12.5e14.2 cm for tiger); pad size width (<6.5 cm for leopard, 8.8e10 cm for tiger), adult stride length (~average of 90 cm in leopard, > 100 cm in tiger), scrapes (<25 cm long and <15 cm wide for leopard, > 35 cm long and >19 cm wide for tiger), and scats (2e4 cm for leopard, > 11 cm in diameter at widest part of the scat for tiger) (McDougal, 1999;Singh et al., 2014;Rostro-García et al., 2016;Simcharoen et al., 2018;Lamichhane et al., 2019). Tiger scat usually is less coiled than leopard scat with a relatively long distance between successive constrictions within a single piece of scat (Wang and Macdonald, 2009). A team of four skilled forest guards and wildlife technicians surveyed each grid cell searching for leopard and tiger signs. There were approximately 100 individuals (~25 wildlife technicians and 75 national park staff, division forest staff, and trained volunteers) involved in surveying after completing training workshops on field survey methods. To avoid biases arising from misclassifications or sign decay (Miller et al., 2011), the team recorded only fresh signs identified with high confidence. To optimize detection, we deviated from search paths (up to 2 m in the case of tracks, to 25 m in case of riverbanks) looking for signs but ensuring uniformity in spatial coverage within each grid cell.

Selection of covariates at the grid cell level
We collected habitat covariates at two spatial scales (landscape and field-based) for predicting occupancy. Most wildlife populations, including leopards, are mostly confined to forested habitat interspersed in the matrix of the human-dominated landscape (Harihar and Pandav, 2012). Available habitat, mostly forest intermixed with grassland, has been defined as potential habitat for leopards in this study and thus governs their presence across the landscape. We extracted landscape-level covariates for each grid cell using a data source downloaded from a geographic information system (GIS) public domain (Table  S1). Covariates included: habitat available (HAB), the extent of deforestation (DEF), and distance to the nearest village (DNV). We calculated habitat availability (HAB) by superimposing the 96-grid cells over the land use and land cover data (Joshi et al., 2003;Reddy et al., 2018) and measured the proportion of habitat within each grid cell divided by total potential habitat (such as sal dominant mixed deciduous forest, grassland, water bodies, etc.). Thus, our occupancy estimates apply only to potential habitat. Within a landscape, high dependency on forest and forest products, encroachments, overgrazing, resettlement, and road building are drivers of deforestation (Kissinger et al., 2012). We measured the extent of deforestation in the landscape. Here, the extent of deforestation is defined as a change in forest canopy cover (in km 2 ) within grid cells between the first (2009) and last (2013) surveys (Thapa et al., 2018). We also conducted forest canopy cover change analysis across the time interval (2009e2013) and extracted this metric for each grid cell. We used 30 m spatial resolution Landsat satellite 7 imagery between September and February for the two-time frames (2009 & 2013). Forest cover classification was analyzed in ERDAS IMAGINE 2015. Post classification and the change analysis were done in Arc GIS (Ver. 10.2). We used standard procedures for forest cover mapping analysis vis-a-vis image acquisition, image stacking, image registration/rectification and masking, supervised classification, and ground verification/truthing followed by an accuracy assessment (Joshi et al., 2003;Reddy et al., 2018).
We used DNV as a surrogate measure of disturbance at the landscape level, which is an easily quantifiable measure that correlates with human activity (Gray and Phan, 2011). We collated settlement point data (~8965 data points) published by the Survey Department (SD, 2001) to compute the distance to the nearest village. The Survey Department spatially defined each settlement point as clusters of households (cluster size undefined) spread across the landscape. We calculated the centroid for each grid cell and computed the DNV using the nearest feature vector-based analytical extension available in ArcView 3.2. We used road abundance (linear km) within each grid cell as a covariate that may affect leopard detection and occupancy across the landscape (Long et al., 2011). We collated road abundance vector data (~9419 km) published by the Survey Department (SD, 2001). From largest to smallest linear road structures, variables included: east-west highway (EWH, paved road), district roads (DR, paved, and/or dirt), and cart roads (CR, dirt road). Details of each type of road network are listed in the supporting information (Table S1). The road (vector) data were superimposed on grid cells and the total length of each road type within each grid cell was summed to obtain road abundance (in km) for each type of road. We expected leopard occupancy to be positively influenced by habitat availability and distance away from nearest villages, and negatively influenced by the extent of deforestation and road abundance (except for cart roads, which could be positive). Comparatively, it is often easier to detect carnivore footprints on cart roads due to soft or muddy substrate conditions, and carnivores readily prefer to use these lowtraffic roads.

Selection of ground-based covariates
At the field level, while walking on trails (1e2 m width) or fire lines (3e4 m in width), the team searched for, and tallied, all signs along the 1 km search paths for each ground level covariate. For each grid cell, the tallies were summed from all 1 km spatial replicates and divided by the total number of spatial replicates.
Poaching, habitat loss, and forest fires are identified as major threats to species in the landscape (MoFSC, 2016). Furthermore, overgrazing by livestock and high resource extraction (for fuelwood, fodder, and non-timber forest products) are threats to habitat at a fine-grain scale. We expected human activities to negatively influence leopards, potentially forcing them to live in sub-optimal habitat with humans, especially outside protected areas (Harihar and Pandav, 2012). We recorded human disturbance signs by counting signs as evidence of livestock presence (sightings, dung, scratches, and tracks). Livestock presence can be considered a surrogate for the intensity of negative human impacts . We also recorded poaching signs (gunshots heard, gun shells found, hunters seen, snares, axes, bottles of poison found, and wildlife carcasses), and evidence of impact on vegetation (lopping and cutting of branches, shrubs, clear felling, presence of logs and/ or carts) along 100 m sections. Fire, as a management tool, is done on a controlled basis and thus it is mostly localized in nature. Wildfires often start in the early winter season (January and February) and continue until late June. We only recorded a fresh evidence of fire that was less than week old as noted by bare ground, smoking embers, and no/little sprouting vegetation. Forest fires are often too large and complicated to quantify in the field. Thus, we considered fire as a binary variable such that presence of forest fires was coded as "1" (present) or "0" for (absent) along every 100 m.
Carnivore densities are known to be a function of prey biomass and density . Thus, local leopard presence should be governed by the availability of wild prey. Wild prey (PREY) presence was based on counting unique signs of all potential prey species combined (sambar, swamp deer, hog deer, chital, barking deer, wild pig) such as dung, tracks, animal sightings, and calls. These allowed us to construct two indices characterizing disturbance (DIST) and wild prey (PREY) that might affect leopard occupancy and detection. Often weighted indices are used based on the perceptions of stakeholders to combine multivariate data into a single index (Salzman, 2003). One such index for quantifying human disturbance was developed by Barber-Meyer et al. (2013). We used this index but note that better methods may be developed in the future that do not rely solely on expert opinion. Nonetheless, we followed Barber-Meyer et al. (2013) in developing a weighted composite index of disturbance (DIST): DIST ¼ 0.35*evidence of poaching þ 0.2*livestock presence þ 0.25* impact on vegetation þ 0.2*fire presence. All the tallied counts (summarized at 1 km) were multiplied by the assigned weights in the equation above. The wild prey index (PREY) was calculated as the sum of prey detections divided by total sampling effort within each grid cell and is expected to positively influence occupancy probability. Also, we modeled survey effort (EFF) as a covariate that might positively influence leopard detection probability, considering that effort varied from cell to cell due to logistical constraints (Harihar and Pandav, 2012) (Table S1). We also predicted a negative influence of road abundance based on road type (EWH and DR) and DIST on leopard detection probability across the landscape. While leopards and tigers are sympatric, there is evidence of spatial and diet partitioning (Seidensticker, 1976;Odden et al., 2010) and we expected to find evidence of spatial segregation at the micro-scale. See supporting information (Table S1) for a complete list of apriori hypotheses.

Data analysis
We used the standard framework (MacKenzie et al., 2002) to model leopard occupancy across the landscape maximizing the likelihood of observing the detection histories at the sites. We used the single species, multi-season model with correlated replicate surveys in Program PRESENCE version 12.2 that explicitly takes into account spatial autocorrelation in detection across the 1 km search paths within each grid cell (Hines et al., 2014). We defined two seasons with the first being 2009 and the second 2013. We used Hines et al. (2014) model to estimate the initial probability of occupancy for the year 2009 and estimated derived probability of occupancy for the year 2013.
All covariates were screened for correlation (Supporting Information, Table S2) and highly significantly correlated variables (i.e. r ! | 0.77|) were either removed or not used in combination within the same model. All covariates used in modeling were normalized using the z-transformation (Sunarto et al., 2012). We tested for autocorrelation by including the spatial dependence parameter (i.e. the autocorrelation parameter at the transect level, q 0 ), which represents the probability that the species is present locally given the species was not present in the previous transect (spatial replicate); while q 1 represents the probability that a species is present locally, given it was present in the previous spatial replicate. We used the constant model for q 0 and q 1 in all analyses. We then used a two-stage approach to model the parameter of interest at the grid cell level . First, we modeled the influence of each of eight covariates (field-DEF, DIST, EFF; landscape-DR, EWH, CR, DNV, HAB) on leopard detection probability, while using the global model (the most parameterized model that included all the landscape and field covariates) for occupancy. After finding the best model for detection probability, we modeled the influence of covariates (landscape-DR, EWH, CR, DNV, HAB; field-PREY, DEF, DIST) on occupancy (Ѱ). We modeled covariates in a stepwise fashion such that if the covariate improved model fit, it was retained in the model and combined with other covariates in multivariate models that we deemed important from our a priori model building. We only used combinations of covariates as additive effects in the models. We eliminated models from the candidate set that did not converge.
We ranked all models using Akaike's Information Criterion (AIC) and chose the best model based on the lowest AIC scores.
We considered all models with DAIC <2 as competing models (Burnham and Anderson, 2002). We used inferences from the 96 surveyed grid cells in winter of 2009 and 2013 to estimate leopard occupancy based on landscape and field level covariates. We modeled site-specific probabilities of leopard occupancy as linear functions of covariates (environmental, habitat, and anthropogenic infrastructures) using the logit link function (MacKenzie et al., 2018). The value of untransformed coefficients (i.e. betas, b), reflects the magnitude and direction (sign, þ/À) of their influences on probabilities of detection and occupancy.
We considered covariates as important and supported if the 95% confidence limits on the b estimates did not include zero (Dupont, 2002). We reported the b estimates for the best model and the univariate models. We computed the model-averaged estimates of grid cell-specific occupancy ( b J) by considering all plausible alternative models by averaging coefficients from competing models based on their AIC model weights. We reported the final occupancy estimates from both the null (constant) model and the top model estimates. To estimate overall leopard occupancy ( c J) in the TAL, we weighed the grid-specific occupancy estimates by potential habitat within each 225 km 2 grid cell . We prepared predictive maps of occupancy for leopards based on inferences made from the model-averaged estimates in ArcGIS 10.1. We also used two species, co-occurrence modeling, accounting for imperfect detection, to analyze species co-occurrence patterns and estimate the species interaction factor, denoted as SIF or f (MacKenzie et al., 2004). We used the 2013 data set for modeling the interactions between tigers (dominant) and leopards (sub-dominant) across the landscape. In PRESENCE version 12.2, we used the psiBa/rBa parameterization, which has increased stability and ability to incorporate covariates in comparison to the phi/delta parameterization (Richmond et al., 2010). Following MacKenzie et al. (2018), we developed two models (one with f estimated and one with f set ¼ 1 denoting independence) and formally compared fit based on AIC values. We considered models competing if DAIC <2 (Burnham and Anderson, 2002). We calculated the species interaction factor (SIF, f), from the top model parameters. A SIF >1 (with CIs that do not overlap 1) indicates co-occurrence more often than expected by chance, while SIF < 1 indicates co-occurrence less often than expected by chance. A SIF ¼ 1 indicates independent occurrence, which is the case when psiBA ¼ psiBa, meaning that species B has the same occupancy whether species A is present (capital A) or absent (little a). We included the covariates deemed important in the single-species, multi-season models (two ecological: habitat available (HAB) and wild prey (PREY); one anthropogenic: disturbance (DIST)) to model the factors affecting the species interaction between leopards and tigers across the landscape. In this way, we can use model rankings to determine whether occupancy was better explained by habitat alone, species interactions alone, or a combination of both. We compared the relationship between prey and habitat on species interaction factors between protected areas and non-protected areas. We conducted linear regression to test the strength (i.e. coefficient of determination (R 2 )) of the influence of habitat and prey on the SIF.

Results
The

Test for spatially correlated occupancy model
The spatial autocorrelation model was the best model as indicated by the lowest AIC values (DAIC ¼ 152, q 0 (SE q 0 ) ¼ 0.13 (0.01),q 1 (SE q 1 ) ¼ 0.73 (0.06), Table 1) and confirmed lack of independence in detections across spatial replicates compared to the standard two-season model. Therefore, we included spatial autocorrelation in all subsequent models.

Effects of covariates on detection and occupancy
We used eight covariates on occupancy probability, Psi (J), and nine covariates on detection probability, p. No continuous variables were correlated, and therefore, all were retained in analyses (Table S2). We compared 28 plausible a priori alternative models (11 for detection and 17 for occupancy), which described expected combinations of the covariates (following Sunarto et al. (2012)) influencing leopard occupancy and detection probabilities. We found the model containing the additive effect of DEF (leopard detection decreased with increase in the extent of deforestation) and district road (DR, leopard detection decreased with increase in abundance of district roads) was the top detection model (w ¼ 0.80, Table 2). The b estimate coefficients for the covariates influencing detection probabilities (Table 3) had varying degrees of influence and confidence intervals (CIs) did not overlap zero indicating support for effects of DEF and DR (Fig. 2), supporting our predictions. There was no other competing model but the next best model (DAIC ¼ 3.52) included DIST such that increasing disturbance led to lower detection probability (w ¼ 0.14, Table 2). Distance to nearest village/settlement did not influence detection, contrary to our predictions. Thereafter, we only modeled the covariates influencing J using the model structure from the top detection models (DEF þ DR) fixed and compared them via AIC.
Three of these competing models contained road type (EWH, DR, CR) impacts influencing occupancy probability (cumulative Table 1 Model support for spatial auto correlation in the multi-season occupancy models incorporating the standard occupancy model with no autocorrelation against the model with autocorrelation (including q 0 , q 1 ) for leopard occupancy across the TAL, Nepal in 2009 and 2013. q 0 (SE q 0 ) ¼ 0.13 (0.01) while q 1 (SE q 1 ) ¼ 0.73 (0.06). J ¼ probability of site occupancy at the grid cell level; p ¼ probability of detection; AIC is Akaike's Information Criterion, DAIC is the difference in AIC value between the top model and the focal model in the set, K is the number of model parameters, and model likelihood is À2 logarithm of the likelihood function evaluated at the maximum. q 0 ¼ spatial dependence parameter representing the probability that the species is present locally, given the species was not present in the previous spatial replicate; q 1 ¼ spatial dependence parameter representing the probability that a species is present locally, given it was present at the previous spatial replicate. g is the probability that the site is occupied in 2009, given that it was unoccupied in 2013. ε is the probability that the site is unoccupied in 2009, given that it was occupied in 2013. pi is the probability that initial replicate is preceded by an occupied site; (.) indicates these parameters are held constant.

Table 2
Effects of covariates on detection probability (p), while using the global model for leopard occupancy using multi-season occupancy models across the TAL, J ¼ probability of site occupancy/habitat use at the grid cell level; p ¼ probability of detection; AIC is Akaike's Information Criterion, DAIC is the difference in AIC value between the top model and the focal model in the set, K is the number of model parameters, and model likelihood is À2 logarithm of the likelihood function evaluated at the maximum. q 0 ¼ spatial dependence parameter representing the probability that the species is present locally, given the species was not present in the previous spatial replicate; q 1 ¼ spatial dependence parameter representing the probability that a species is present locally, given it was present at the previous spatial replicate. g is the probability that the site is occupied in 2013, given that it was unoccupied in 2009. ε is the probability that the site is unoccupied in 2013, given that it was occupied in 2009. pi is the probability that initial replicate is preceded by an occupied site. Covariates considered: PREY: wild prey index, DR: district road, DEF: extent of deforestation, DIST: anthropogenic disturbance index, EFF: realized survey effort, EWH: east west highway, DNV: distance to nearest village, CR: cart road, HAB: total habitat available in each of grid cell. (.) denotes parameter held constant. In all models, occupancy, "J", was modeled as Global:  Table 3 The b estimates and standard errors (in parentheses) from the logit link function based on the best, competing, and the eight univariate multi-season occupancy models for leopard detection probability (p) in TAL, Nepal in 2009 and 2013. b estimate is presented for each of the univariate models where each model represents one covariate.  "C" represents mean; "C" represent upper confidence limit; ":" represent lower confidence limit; DEF: extent of deforestation; DR: district road.

Table 4
Effects of covariates on leopard occupancy (J) while using the top two detection models across the TAL, Nepal from multi-season occupancy models in 2009 and 2013. Also included is the null (constant) model for comparison. J ¼ probability of site occupancy/habitat use at the grid cell level; p ¼ probability of detection; AIC is Akaike's Information Criterion, DAIC is the difference in AIC value of the focal model and the best AIC model in the set, K is the number of model parameters; model likelihood À2 logarithm of the likelihood function evaluated at the maximum. q 0 ¼ spatial dependence parameter representing the probability that the species is present locally, given the species was not present in the previous spatial replicate; q 1 ¼ spatial dependence parameter representing the probability that a species is present locally, given it was present at the previous spatial replicate. g is the probability that the site is occupied in 2009, given that it was unoccupied in 2013. ε is the probability that the site is unoccupied in 2009, given that it was occupied in 2013. pi is the probability that initial replicate is preceded by an occupied site. Covariates considered: DEF: extent of deforestation; DIST: anthropogenic disturbance index; DNV: distance to nearest village; HAB: total habitat available in each of grid; PREY: wild prey index; EWH: east west highway; DR: district road; CR: cart road; (.) denotes parameter held constant. In all models, the probability of detection (top model from Table 2) was modeled as "p 1 " (PREY þ DR)". 'þ' denotes covariates modeled additively. Post Hoc Analysis: p 2 "PREY þ DEF" as second-best detection model (Table 2) was added to evaluate influence on occupancy.

DAIC
AIC weight ¼ 0.17; DAIC<2; Table 4). Except for cart road (CR), which was positive, the larger east-west road (EWH) and district road (DR) abundances had negative effects on leopard occupancy following our a priori predictions, however, CIs overlapped zero indicating inconclusive effects (Table 5). Deforestation (DEF) and disturbance (DIST) negatively affected leopard occupancy while increasing distance to the nearest village (DNV) increased leopard occupancy in concordance with our a priori predictions, but again, CIs overlapped zero (Table 5). As expected, habitat patch size and especially prey index, positively influenced leopard occupancy across the landscape, but CIs overlapped zero ( Fig. 3; Table 5).  The b estimates and standard errors (in parentheses) from the logit link function based on the best and the univariate, multi-season occupancy models for leopard occupancy (J) in TAL, Nepal, in 2009 and b estimate is presented for each of the univariate models where each model represents one covariate. DEF: extent of deforestation; DIST: anthropogenic disturbance index; DNV: distance to nearest village; HAB: total habitat available in each of grid, Prey: wild prey index, EWH: east west highway; DR: district road; CR: cart road; SE: standard error. Bold indicates strong or robust impact, that is 95% confidence intervals as defined by b b±1:96 Â SE not overlapping at 0; Italics indicate opposite from a priori prediction. Fig. 3. Relationships between continuous covariates and leopard occupancy from the competing univariate, multi-season models in TAL, Nepal, in 2009 and 2013. "C" represents mean; "C" represent upper confidence limit; ":" represent lower confidence limit; DEF: extent of deforestation; HAB: habitat available; EWH: east west highway; DNV: distance to nearest village; PREY: wild prey index; DR: district road; DIST: disturbance; CR: cart road.  Two species co-occurrence model results including test for species interaction (in italics; DAIC > 15.09) between tigers and leopards in the TAL, Nepal from 2013. total 13,845 km 2 of potential available habitat in the TAL was 12,184 km 2 (SE 692 km 2 ) for 2009 and dropped to 8168 km 2 (SE 831 km 2 ) in 2013. We estimated grid cell-specific b J and variation across geographical space in the TAL, Nepal (Fig. 4).

Leopard-tiger interaction at the landscape scale
Based on the delta AIC between models that estimated the SIF and those with SIF fixed at 1.0 (independence), we found strong statistical support for an interaction between tigers and leopards (DAIC > 15.09, Table 6). Based on estimated SIF of 1.49 (SE 0.13), we found spatial co-occurrence between tigers and leopards rather than avoidance at the landscape level. We also found SIF to be meaningfully different from 1.0 (i.e., the 95% confidence intervals did not include 1.0) across all grid cells. Model-averaged occupancy of leopards when tigers were present was higher (0.85 ± 0.04) than when tigers were absent (0.39 ± 0.01), as expected under co-occurrence, BA > Ba. However, this relationship was mediated by habitat patch size (Table  6). In grid cells where tigers were detected, leopard occupancy decreased with increasing habitat availability; while in grid Fig. 5. Relationship between habitat and leopard occupancy (with 95% of CI) in the presence of (psiBA) and in the absence of (psiBa) tigers based on top model (top left); relationship between species interaction factor (with 95% of CI) and habitat availability (top right). Species interaction factor of 1.0, indicating independence, is denoted by the dotted line. Influence of habitat and prey on species interaction factor between grids in protected areas (middle and bottom left) and area outside of protected areas (middle and bottom right) along with coefficient of determination (R 2 ) for regression analysis. cells where tigers were not detected, leopard occupancy increased with increase habitat availability. Interestingly, while the SIF was greater than 1.0, it decreased as habitat availability increased, approaching independence at very large habitat patches (Fig. 5).
Within protected areas, we found that tigers and leopards co-occurred in small habitat patches, but as patch size increased to~>75 km 2 , they occurred independently. We also found leopards and tigers strongly co-occurred at low prey index values within protected areas, likely indicating they are attracted to the same areas when prey is scarce, while they occurred independently at moderate to high prey index values. Outside protected areas, the SIF declined more strongly as patch size increase (as indicated by regression coefficient), but it never reached independence and always indicated co-occurrence between tigers and leopards even in large habitat patches. The prey index values were much lower outside protected areas and, in contrast to protected areas, leopards and tigers always co-occurred regardless of prey index (although the SIF declined slightly), likely because of the overall scarcity of prey (Fig. 5).

Discussion
Species distribution is a fundamentally important parameter in ecology for informing large-scale conservation and management (Jones, 2011;Karanth et al., 2011). The use of multi-season, robust design occupancy models provide an opportunity to empirically identify the distribution of multiple carnivores on the landscape. Our results indicate that a) leopard occupancy declined by 31% between 2009 (Psi: 0.86) and 2013 (Psi: 0.59), b) leopard detection was strongly, negatively influenced by deforestation and district roads, c) anthropogenic factors such as disturbance, road network (particularly the east-west highway), and extent of deforestation (decrease in forest cover) appear to be negatively affecting occupancy across the landscape, but effects were generally weak or inconclusive; d) tigers and leopards tend to co-occur across the landscape but this interaction is mediated by extent of available habitat and prey index, in slightly different ways inside and outside protected areas.
We note that signs (track, scats, scrapes) of young tigers can overlap those of leopards in size, however, young leopardsized tigers were likely rare and primarily limited to sites within protected areas. So mistaken identification of tigers as leopards would be few. Additionally, these errors would be unlikely to alter any of our inferences on occupancy since areas where young tigers would be unlikely to alter our inferences on occupancy since areas where young tigers were likely to occur, were always occupied by leopards. These potential errors in misidentification could inflate our detection probability for leopards, but this error would be consistent error for surveys in both time periods.
There was support for the Hines et al. (2014), spatial autocorrelation models, compared to the standard occupancy models (MacKenzie et al. 2.002) indicating a lack of independence in sign detection on search paths along the 1 km spatial replicates. Thus, we could test sets of covariates, including ecological and anthropogenic variables, while incorporating spatial autocorrelation.
We found that the probability of leopard detection was negatively related to the length of district roads (expressed in km per grid cell) and the extent of deforestation (in km 2 ) in each grid cell across the landscape. Our results also indicated that leopards tend to avoid the larger district roads. Thus, the construction of large rural roads (asphalt or larger dirt roads) including a high density of roads in potential habitat is likely to have negative impacts on leopards. The extent of deforestation was found to have a strong negative influence on leopard detection. Encroachment into the potential forested habitat is observed as one of the major drivers of deforestation (MoFSC, 2016) on the landscape, which may, therefore, affect detection probability.
For occupancy, though all covariates (in competing univariate models) conformed to our a priori expectations, none of the covariates was strongly supported (CIs overlapped 0). This could be due to the generalist nature of the leopard itself, or perhaps the degree of specificity of covariates used in the analyses. As a potential weakness in our design, the sampling grid cells may have been too large for the covariates to be strong predictors of leopard occupancy. A potentially more appropriate design specific to leopard home range size could be~10 km Â 10 km grid cells (or even smaller sizes) with covariates also collected on this smaller scale as done in India (Gubbi et al., 2020). However, the surveys were conducted following the same protocol in 2009 and 2013, and we still documented a dramatic decline in leopard occupancy, despite the large grid cell size. We additionally gained insight on potential drivers of occupancy, even with the potentially low power of co-variates and we capitalized on the opportunity to analyze data for a species of concern sympatric to tigers.
Although all effects were inconclusive, competing models for occupancy included anthropogenic factors (DEF, DNV, DIST; cumulative AIC weight (cw): 26%), ecological covariates (PREY and HAB; cw: 13%), and road abundance (EWH, DR, CT; cw:18%). Though leopards are believed to be able to live close to people, human-caused disturbances and infrastructure development were the strongest predictors causing declines in leopard occupancy on the landscape.  stressed the adoption of strategies by leopards that minimize direct contact with humans in high human-disturbed areas. However, due to the lack of fine-grain, spatio-temporal patterns in our study area, more research is needed to link occupancy with drivers of the decline. Similar results were found by Kshettry et al. (2017) where leopard habitat use was low and negatively influenced by high the density of infrastructure (such as buildings) in the Western Bengal landscape, India.
Large linear infrastructures are a direct threat to the loss of habitat in the landscape (MoFSC, 2016). Among road networks, highways and district roads appear to be affecting (although weakly) leopards at the landscape level. Our results support earlier studies that road networks are impacting wildlife roaming around their edge: Amur tigers in Russian far East (Kerley et al., 2002), large mammals in India (Gubbi et al., 2012), elephants in Malaysia (Wadey et al., 2018), and tigers in their range countries (Carter et al., 2020). The Government of Nepal (GoN) constructed a total of 57,632 km of roads across the country, a 23-fold increase from 2500 km since the late 1950s (Bhandari et al., 2012). The GoN revealed that the fifteenth five-year strategic plan expects to construct more district roads and national highways to increase national accessibility by 2030 (MoPPT, 2017). Among the suite of mitigation measures available, wildlife crossing structures are planned at strategic sites where wildlife crossing is evident (Fig. 6). However, as precautionary measures against loss of biodiversity, development efforts should avoid potential wildlife habitat to the extent possible. Mega-development projects such as the widening of the NH7 highway in India (Dutta et al., 2013) and the East-West highway in Nepal  were directed to strictly follow mitigation measures that facilitate dispersal, mediate geneflow, and enforce measures to avoid vehicular collisions for a suite of mammals. While we only found weak impacts of the road network at the current time, all impacts of these larger road networks (East-West Highway and District Roads) on leopard occupancy were negative while small cart roads (CR) positively influenced occupancy indicating that leopards are not deterred by such small roads and likely use them for travel. Continuing multi-season occupancy analysis in the future may reveal stronger negative effects as more infrastructure is built in Nepal. We may be seeing the beginning of what are likely to be strong negative impacts of large roads on future leopard occupancy.
The large decline in occupancy over time could be attributed to habitat loss and destruction, poaching, and retaliatory killing arising from human-leopard interactions. A recent nationwide forest assessment showed that lowland Terai forests have been depleted due to deforestation at an annual rate at 0.44% (16.5 km 2 /year) between 2001 and 2010, whereas forest cover has slightly increased in mid-hills (DFRS, 2015). As a caution, the leopard was not found in many sites that were occupied in 2009 including some sites within protected areas. Analysis of the records maintained by the Department of National Parks and Wildlife Conservation (DNPWC) highlighted that leopard pelt seizures increased by threefold between 2009 and 2013 (Paudel et al., 2020). Due to improved strengthening of protection units within the protected areas system, poaching of leopards is now minimal but it is a concern outside protected areas. Thapa (2015) reported human-induced mortality (poaching, lethal control, retaliatory killing, and road accidents) at~63% with 45 individuals known to be killed between 2009 and 2013 in Nepal. However, the actual mortality of leopards likely is higher as many cases are not reported. Both mortality and illegal trade records suggest poaching of leopards is high in the region. In India, annually 155 leopards are killed due to poaching alone (WPSI, 2018). This could explain the decrease in leopard occupancy across these five years in the TAL, Nepal.
Leopards were colonizing grid cells in 2013 that were not occupied in 2009, especially in grids cells outside protected areas. This shows leopard may have strongholds in protected areas and maybe attempting dispersal into surrounding areas with suitable habitat. Interestingly, while leopard occupancy declined over the period, tiger occupancy increased (ѱ 2009 -0.37 to ѱ 2013 -0.55) across the landscape (Barber-Meyer et al., 2013;Dhakal et al., 2014;DNPWC & DFSC, 2018). As a caveat, we did not present observed dynamic occupancy from two season only as a cue to observe trend in leopard occupancy. Such analysis would require more seasons of data to assess variation in this dynamic occupancy to make stronger conclusions about the strength in the overall occupancy trend and serves important implications for future research.
Studies within protected areas in Nepal, Chitwan (Seidensticker, 1976), Bardia (Odden et al., 2010), and Shuklaphanta (Lovari et al., 2015), which collectively explain mechanisms of co-occurrence at a fine grain scale, suggest that habitat partitioning occurs through the activity of sympatric carnivores over space, time, and through their food habits. While we could not include temporal partitioning or diet in our analyses, we found that tigers and leopards co-occurred at the landscape level. In contrast, Kafley et al. (2019) found spatial avoidance between the two predators at a fine scale (2 km Â 2 km grid cells) in Chitwan National Park based on a camera trap study. This indicates that mechanisms governing co-occurrence could differ between scales. We found co-occurrence between tigers and leopards was facilitated by habitat availability and wild prey at the landscape level rather than being driven by negative interactions expected via interspecific competition.
Within protected areas, tigers and leopards co-occur in small habitat patches when prey was low, presumably because they are attracted to a scarce resource in a small space. But as patch size or prey index values increase, they occur independently, implying they no longer need to use the same areas to meet their spatial requirements or energy demands. We did not see this outside of protected areas. Instead, tigers and leopards always co-occurred, even at large patch sizes (defined here as maximum of 225 km 2 grid cell), and regardless of prey index. Although never reaching independent occurrence, the SIF strongly declined with patch size, potentially indicating that, outside of protected areas, patch sizes need to be larger than is currently available on the landscape for the two species to occur independently. Co-occurrence between leopards and tigers was especially high in small habitat patches outside protected areas, indicating that these areas could be dangerous for humans and potential hotspots of human-wildlife conflict. Additionally, the prey index outside of protected areas was always lower (by an order of magnitude) than inside protected areas. Thus, it is likely tigers and leopards also co-occur outside protected areas due to the need to find prey that is always scarce.
Our analyses were conducted on a broad scale. We suggest additional monitoring of carnivores at a finer scale through intensive camera trapping, fecal DNA surveys, and radio telemetry studies to provide more refined insight. In future surveys, we plan to analyze scats genetically to estimate landscape-level abundance and density and to explore genetic connectivity across the landscape. The GoN has the initiative to monitor carnivores using camera trapping at fine-scale levels and this offers an opportunity to integrate the occupancy modeling framework with the spatially-explicit, mark-recapture density estimation framework (Royle et al., 2013). We know that leopard density (no of individuals per 100 km 2 ) varies with habitat type: churia (4.00 (Thapa and Kelly, 2017),), bhabar (3.48 ,); lowland floodplain (3.12 (Thapa, 2011),). Yet, it is still a challenge logistically and financially to survey the entire landscape at a fine-scale to estimate leopard density. Also, to continue multi-season occupancy modeling to track distributional changes over time, we suggest that future research utilize recent advances (Dey et al., 2017;Linden et al., 2017) integrating data from multiple sources at multiple scales, such as camera trap data from core areas (Kumar et al., 2019) and landscape-wide occupancy data, as a viable way forward for estimating landscape-scale leopard density.

Conclusion
This study is the first of its kind to analyze leopard occupancy along the spatially replicated transects designed for tigers and leopard interactions with a congener at the landscape level. Tigers and leopards were found to co-occur across the landscape, but the amount of available habitat and prey mediated co-occurrence. Mitigation measures should focus on preventing loss of critical leopard, tiger, and prey habitat through appropriate wildlife-friendly underpasses and avoid such habitat when building infrastructure. The GoN flagship monitoring program should continue to monitor large carnivore using the standard protocol across the Terai Arc Landscape and beyond in mid-hills and high mountains. Repeated surveys through time will provide data on occupancy trends that are useful for protected area managers and natural resource management communities (e.g. community forest user groups) for devising appropriate management strategies on the ground. Leopard conservation has received a lower priority than tigers, but our metrics show a large decline in leopard occupancy over time. Conservation planning should focus on measures to reverse this decline by facilitating human-leopard coexistence to ensure leopard persistence across the landscape.

Author contributions
KT, SM, SRJ, NS, MD conceived the study; KPA, MD, MTK, SRJ, SRB administered the study; KT, SM, SAS, NS, BRL collected the data; KT and MJK analyzed the data and wrote the paper; all authors reviewed the paper.

Funding
This work was funded as a part of support received as part of tiger monitoring project from WWF Network (UK, US and International), US Fish and Wildlife Service, USAID's Hariyo Ban Program through WWF Nepal, The Leonardo DiCaprio Foundation, and Save the Tiger Fund. Virginia Tech for its support to KT during the analysis and manuscript development.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.