Tigers on the edge: mortality and landscape change dominate individual-based spatially-explicit simulations of a small tiger population

Reductions in the tiger’s (Panthera tigris) range in Southeast Asia have been concurrent with large infrastructure expansion and landscape change. Thailand’s Dong Phayayen-Khao Yai Forest Complex (DPKY), a landscape of tiger conservation priority, may be particularly vulnerable to such changes, necessitating investigations into effects on population dynamics. Evaluate relative effects of landscape change scenarios on the probability of tiger persistence in DPKY and sensitivity of predictions to spatially-explicit mortality risk, landscape resistance, and tiger population density. We utilize individual-based, spatially-explicit population modelling to evaluate the trajectory of tiger population dynamics across 11 landscape change scenarios. Concurrently, we evaluate sensitivity of predictions to landscape resistance transformation, maximum population density, and spatially-explicit mortality across 20 generations. Spatially-explicit mortality risk dominated predictions of population persistence, frequently resulting in population declines/extinction. Adjustment of moderate mortality risk to slightly convex and concave forms shifted extinction rates from 46 to 12% and 85%, respectively. Holding mortality constant at moderate levels, strong negative effects were predicted in landscape change scenarios incorporating road expansion (46%-74% extinction) and construction of dams (52%). Strong negative effects of combined development persisted even when habitat restoration measures were applied (96% extinction). Adjusting resistance and maximum population density had marginal effects. The high sensitivity and variability of predictions to spatial patterns of mortality risk suggest a population on a proverbial knife’s edge. Our results underscore the importance of incorporating spatial patterns of mortality risk in population modelling, highlighting their potentially dominating influence on population dynamics and extinction risk.

spatially-explicit mortality risk, landscape resistance, and tiger population density. Methods We utilize individual-based, spatiallyexplicit population modelling to evaluate the trajectory of tiger population dynamics across 11 landscape change scenarios. Concurrently, we evaluate sensitivity of predictions to landscape resistance transformation, maximum population density, and spatiallyexplicit mortality across 20 generations. Results Spatially-explicit mortality risk dominated predictions of population persistence, frequently resulting in population declines/extinction. Adjustment of moderate mortality risk to slightly convex and concave forms shifted extinction rates from 46 to 12% and 85%, respectively. Holding mortality constant at moderate levels, strong negative effects were predicted in landscape change scenarios incorporating road expansion (46%-74% extinction) and construction of dams (52%). Strong negative effects of combined development persisted even when habitat restoration measures were applied (96% extinction). Adjusting resistance and maximum population density had marginal effects. Conclusions The high sensitivity and variability of predictions to spatial patterns of mortality risk suggest a population on a proverbial knife's edge. Our results underscore the importance of incorporating spatial patterns of mortality risk in population modelling, highlighting their potentially dominating influence on population dynamics and extinction risk.

Introduction
Tigers (Panthera tigris), which historically occurred across Asia, have suffered large contractions in range and population over the past century (Nowell and Jackson 1996;Goodrich et al. 2015). As human populations continue to expand, refugia for tigers are becoming increasingly scarce and fragmented. This is compounded by poaching and other threats, jeopardizing long-term population persistence (Wikramanayake et al. 2011;Goodrich et al. 2015;Wolf and Ripple 2016). This is particularly evident in Southeast Asia. Severe reductions in the tiger's range in the region (Sanderson et al. 2006;Lynam and Nowell 2011) have occurred at a time where national and transnational economic strategies have necessitated monumental investments in infrastructure and transformative land use change (Quintero et al. 2010;Wongwuttiwat and Lawanna 2018;Hughes 2019;Kaszta et al. 2020). With tigers in Southeast Asia seemingly at both figurative and, often literal, crossroads, understanding the effects of these changes on tiger population dynamics is crucial for developing effective conservation and management strategies.
Given anthropogenically-driven population and range declines over the past century, population dynamics and persistence of tigers currently is less a reflection of the species' biological or adaptive constraints than an indication of spatial patterns of human presence and activity (Sunquist et al. 1999;Walston et al. 2010;Goodrich et al. 2015). Observed demographic patterns governing species persistence are inexorably linked with spatially-dynamic processes such as large scale-changes to landscapes (Lande 1988;Kareiva and Wennergren 1995;Hanski 2011), which is of particular relevance for wideranging species such as the tiger. Understanding these links between species persistence and heterogeneous landscapes can be aided by models which incorporate a spatially-explicit framework (Kareiva 1990;Dunning et al. 1992;DeAngelis and Yurek 2017). Spatially-explicit models can enable evaluation of the effect of changes in the landscape on population dynamics (Dunning Jr et al. 1995;Wiegand et al. 2005), elucidate interactions between landscape structure, physiognomy, and species traits (Cantrell and Cosner 1993;Watkins et al. 2015), and identify potential extinction thresholds (Kareiva and Wennergren 1995;Bartlett et al. 2016).
Several recent studies have conducted spatiallyexplicit simulations of a range of realistic development scenarios and conservation interventions, assessing their effect on vulnerable species (Cushman et al. 2016;Thatte et al. 2018;Macdonald et al. 2018;Kaszta et al. 2019Kaszta et al. , 2020. Importantly, when spatiallyheterogeneous mortality risk is included in these simulations, it can exert a disproportionate influence on population persistence (Kramer-Schadt et al. 2004;Cushman 2006;Kaszta et al. 2019Kaszta et al. , 2020. Despite its potential effect, studies may not explicitly incorporate mortality in modelling relationships between species populations and changes to landscapes (Thatte et al. 2018;Kaszta et al. 2019). It is crucial to investigate not only the effect of changes in landscape configuration on population connectivity and habitat, but also how associated changes in the spatial pattern of mortality risk affect the population and its probability of extinction. Predictions of population dynamics across landscapes may also be sensitive to population size or density (Cushman et al. 2013;Ash et al. 2020a). Population density may influence movement of species across landscapes through availability and selection of habitat (Pulliam and Danielson 1991) though relationships between population density and dispersal in carnivores may be complex (Støen et al. 2006;Joshi et al. 2013). Further, the degree to which higher relative population densities may mitigate effects of mortality in large carnivores is subject to speculation (Karanth and Stith 1999;Chapron et al 2008). Rapid land use change and economic development across the tiger's range also have implications for landscape resistance to movement, which may have drastic and complex effects for highly vagile species (Carr and Fahrig 2001;Carr et al. 2002;Cushman 2006;Cushman et al. 2010). Importantly, predictions may be sensitive to both the means by which factors in population models are parameterized and the strength of their interactions (Dunning Jr et al. 1995;McCarthy et al. 1995).
Thailand's Dong Phayayen-Khao Yai Forest Complex (DPKY) is both an illustrative proxy for Southeast Asia's tigers and is a dynamic landscape for which spatially-explicit simulations may be valuable.

3
Vol.: (0123456789) The UNESCO World Heritage Site supports one of the region's remaining tiger breeding populations (Ash et al. 2020c), though the effects of human pressure on the landscape has been the subject of concern (IUCN World Heritage Outlook 2020). The landscape contains large extents of forest cover of varying contiguity and protection status, with potential opportunities to expand tiger habitat. However, the complex remains entrenched within an increasingly human-dominated landscape at the nexus of important economic corridors (Wongwuttiwat and Lawanna 2018). Roads in and around the complex have and may undergo further expansion and increased use while the placement of wildlife crossing structures along highways could potentially mitigate road effects on population fragmentation (IUCN World Heritage Outlook 2020; Paansri et al. 2021). Lastly, growing concerns over local water insecurity could increase pressure to construct new dams/reservoirs, potentially inundating parts of the complex (Marks 2011;IUCN World Heritage Outlook 2020;Manorom 2020).
Given the importance of DPKY within regional tiger conservation strategies, investigating the potential effect of future development and conservation initiatives is of critical importance. In this study, we utilize an individual-based, spatially-explicit population modelling approach to evaluate the relative effects of potential landscape change scenarios on the probability of tiger persistence in DPKY. Concurrently, we evaluate the relative effect of spatially-differential mortality risk, landscape resistance transformation, and maximum potential population density, as well as the sensitivity of predictions to these factors. Through this study, we aim to generate comprehensive insight for the development of conservation and management strategies in DPKY while highlighting key considerations for spatially-explicit population modelling for other threatened species.

Study site
Eastern Thailand's Dong Phayayen-Khao Yai Forest Complex (DPKY) consists of five protected areas across 6,155km 2 -Khao Yai National Park, Thap Lan National Park, Pang Sida National Park, Ta Phraya National Park and Dong Yai Wildlife Sanctuary -and supports a small, isolated breeding population of Indochinese tigers (Panthera tigris corbetti). The complex is primarily surrounded by a human-dominated landscape of urban areas, villages, agriculture, roads, and other infrastructure. While most of the landscape outside protected areas in DPKY has been converted by human activities, limited patches of forest occur outside the complex adjacent to existing protected areas. The complex is divided by major highways, including routes 304 and 348, though the former has been the site of a recently completed series of wildlife crossing structures, potentially mitigating its fragmentary effects. Restricted or managed roads include a private road (Route 3446) running through Ta Phraya National Park, and an unpaved former public road (Route 3462) which runs through Pang Sida and Thap Lan national parks, managed by park authorities. Outside the complex, minor roads connecting villages and urban areas are abundant. DPKY is the source of major watersheds and a number of man-made dams and reservoirs which primarily occur at the boundaries of, but, in some cases within, protected areas. These are typically maintained for the purposes of local water security and use, such as for agricultural irrigation. The complex is also relatively close to developed urban areas of varying size (including Saraburi [~ 28 km] and Nakhon Ratchasima [~ 40 km]) and, due to its proximity to the Cambodian border, is close to both formal and informal border crossings (Fig. 1).

Population modelling approach
We evaluate tiger population trajectories using CDPOP (Cost Distance POPulations; Landguth and Cushman 2010) across combinations of four key factors: (1) landscape change scenario (SCN); (2) resistance surface transformation (RST); (3) maximum population density (DEN); and (4) mortality function (MRT; see sections below). We vary each factor by several levels, described below, to evaluate the relative effect of these factors on the resulting simulated population (N). A workflow diagram summarizing this process can be found in Fig. 2.

Base resistance definition
In order to test the effect of landscape change and define spatial patterns of the study landscape with which individuals in simulations would interact, we developed a suite of resistance surfaces. Resistance surfaces are commonly used in modelling species movement across dynamic landscapes (Zeller et al. 2012) and reflect step-wise cost to individual movement, with high pixel values representing high resistance to movement and low values imparting little resistance. In our study, resistance surfaces were developed at a 250 m resolution, using ArcGIS 10.3.1 (Environmental Systems Research Incorporated, ESRI, Redlands, CA, USA, 2011) as well as the raster package in R (R Development Core Team 2019; Hijmans 2020).
We derived current landscape resistance from an expert-based resistance surface developed for DPKY by Ash et al. (2020a) who identified the resistance model used in their paper as the most representative of likely patterns of resistance among four candidate models. This expert-based surface was based on 2018 SERVIR-Mekong Regional Land Cover Monitoring System (RLCMS; SERVIR-Mekong 2018) land cover data, classified into dense forest (resistance value of 1), scrub forest (resistance of 20), agriculture/village matrix (resistance of 50), reservoirs/surface water (resistance of 80), and urban areas (resistance of 100). Minor (resistance of 30) and major roads (resistance of 100; OpenStreetMap 2019) were also included. We utilized this base resistance surface as a foundation for the development of further landscape change scenarios.
In order to account for a broad range of potential changes to the landscape, in addition to the base resistance surface described above (herein described as B0), we developed additional resistance surfaces based on a range of landscape change scenarios (SCN; Table 1). In lieu of specific development and landscape change plans, we explored a range of potential and hypothetical changes to this landscape. These included five development scenarios (D1-D5) and two restoration scenarios (R1 and R2). Further, in addition to developing landscape change scenarios independently, we developed three additional resistance surfaces to test the collective effects of changes (D0, R0, DR0). These landscape changes were also included in sensitivity analysis of resistance surfaces. While the extent to which such changes will be realized is uncertain, proactive assessments of potential impacts can be valuable in understanding likely effects prior to changes taking place (Cushman et al. 2016). This is particularly  important given such changes, particularly development, can occur rapidly and without warning.

Development scenarios
We tested potential effects of increased infrastructure development and human activity (e.g., agricultural or urban expansion) were tested via a range of development scenarios. In order to investigate the potential importance of forest cover in and around DPKY for tiger movement and persistence, we developed a resistance surface (D1) in which all existing forest cover outside protected areas was converted to comparatively high-resistance agriculture/village matrix (resistance of 50; sensu Cushman et al. (2016).
In scenario D2, we evaluated potentially increased effects of roads in and around DPKY by doubling the original resistance of both minor (new resistance of 60) and major (new resistance of 200) roads. Further, in scenario D3, we tested the impact of a potential re-opening of decommissioned Route 3462 (IUCN World Heritage Outlook 2020) by All scenarios combined All development and restoration scenarios described above digitizing the former public road and with a resistance value of 100. Scenario D4 included 7 proposed dams built at major rivers throughout the DPKY forest complex which could potentially inundate parts of protected areas (IUCN World Heritage Outlook 2020). Dams and subsequent inundation zones were simulated via ArcGIS 10.3.1, with dam heights corresponding to the recently completed Huay Samong Dam (Naruebodindrachinta reservoir; ~ 33 m). At each approximate dam location, we extracted a section of digital elevation surface, from a minimum elevation at the proposed location up to an additional 29 m, simulating a potential inundation area from a dam height of 33 m with a supply level 4 m below capacity. The resulting surface was treated as a potential inundation zone and assigned a resistance of 80, matching existing reservoirs in the base resistance surface.
In scenario D5, we aimed to simulate hypothetical expansion of urban centres and border areas as a result of economic development and increased trade with neighbouring Cambodia (Kobayashi et al. 2017;Wongwuttiwat and Lawanna 2018). We simulated growth by applying the 'Expand' function in ArcGIS. Each urban pixel was expanded by 1 into neighbouring cells, roughly doubling total urban extent through expansion of existing urban areas. We also simulated the development of hypothetical border crossings at Ta Phraya NP in which urban pixels were generated randomly at two candidate border locations. Pixels were generated as a function of distance from border road with the number of pixels generated capped at the number present at the nearest major border crossing (Aranyaprathet-Poipet). We simulated associated increases in road traffic at these border crossings were simulated by increasing resistance values of private road Route 3446 and other connecting minor roads from 30 (B0) to 100. Further, to simulate the potential diffusive effect of increased urbanization, we applied a 5,000 m point density buffer from all urban pixels (sensu Elliot et al. (2014) declining linearly from a max resistance of 50.

Restoration scenarios
In addition to development scenarios, we also developed two restoration scenarios to evaluate the potential effects of forest reforestation, expansion of protected areas (scenario R1), and construction of wildlife crossing structures on key roads (scenario R2).
In scenario R1, we identified several areas for potential reforestation, primarily adjacent to existing protected areas, were identified using satellite imagery (Google, US Dept of State Geographer, Image Landsat/Copernicus) and forest quality records generated during camera-trap field surveys (Ash et al. 2020c). Resistance values in these areas were reclassified with a resistance value of 1. In addition, to examine the potential effect of increased protected area coverage, where existing or reforested areas occurred outside and adjacent to PAs, we converted forest coverage to polygons and merged them with protected area shapefiles (UNEP-WCMC and IUCN 2019) as hypothetical protected areas extensions. This resulted in an approximate increase of protected area coverage by ~ 14% (+ 887km 2 ) and forest cover by ~ 212km 2 .
In scenario R2, we tested the addition of wildlife crossing structures along two major roads in the DPKY forest complex, Route 304 and Route 348 (Ash et al. 2020a). This includes a highway overpass (~ 570 m) and two tunnels in southern DPKY (~ 210 m and ~ 150 m) and two overpasses in northern DPKY (~ 260 m and ~ 180 m). Though not currently planned, we also included a hypothetical crossing structure located along Route 348 where the highway crosses a forested section of Ta Phraya National Park (one overpass of ~ 570 m). At these locations, we reclassified the resistance of the area occupied by crossing structures as 20 (corresponding to the resistance value of scrub forest).

Combined scenarios
In addition to developing landscape change scenarios independently, we developed three additional resistance surfaces to test the potential additive and interactive effects of these changes. These included: (1) Scenario D0, which reflected full development by incorporating all landscape changes of development scenarios (D1-D5); (2) Scenario R0, which included all restoration scenarios (R1 and R2); and (3) Scenario DR0, incorporating all development and restoration changes. Landscape metrics quantifying the degree of landscape change under these scenarios are presented in Table 2.

Resistance transformation (RST)
Parameterization of resistance surfaces in population connectivity and related studies is ideally informed by movement, genetics, and other empirical data (Zeller et al. 2012). Given the lack of availability of these data, inherent uncertainty in developing resistance surfaces, and accounting for potential non-linear relationships with resistance, we modified resistance for each landscape change scenario via three power function resistance transformations (RST; sensu Reddy et al. 2017). We transformed base resistance (RST10) to the 0.5 (RST05) and 2.0 power (RST20). We then rescaled the result to original minimum and maximum values and resampled them to 250 m resolution. These power functions adjust the intensity and form of resistance from linear (RST10) to concave (RST05) or convex (RST20) functional forms.

Population density (DEN)
Population density (DEN) in our study was defined by a grid of occupiable locations across our study extent. This grid was generated from a cost-distance kernel of 45,000 units from DPKY's contiguous forest boundary which was resampled to 10 km resolution and converted to points. This effectively distributed potentially occupiable tiger home ranges at a theoretical population density of 1/100 km 2 , occurring in the upper range of population density estimates for the DPKY forest complex (Ash et al. 2020b). In order to evaluate the effect of tiger population density on predicted population trajectory, we also tested two additional maximum tiger densities. This includes points generated at a 7.5 km resolution, corresponding to a tiger population density of approximately 1.78/100km 2 , within a range of estimates reported in western Thailand (Duangchantrasiri et al. 2016). We also generated points at 6 km resolution, corresponding to a relatively high theoretical density of 2.78/100km 2 . A total of 104 points were generated at 10 km, 187 points at 7.5 km, and 291 points at 6 km, defining area-specific habitat carrying capacity (K).

Initial Simulations in CDPOP
We ran simulations in CDPOP, which uses an individual-based, spatially-explicit framework to model mating, dispersal and mortality through space and time. Individual movement is defined by sex-dependent movement parameters, differentiated for mating and dispersal behaviour. CDPOP simulates mating among individuals based on probabilistic selection of mating pairs based on the cost distance between their locations across the resistance surface. Dispersal of offspring is also probabilistic based on cost distance, governed by the specified dispersal function and maximum cost threshold. These and other simulation Table 2 Landscape metrics of dense forest (resistance value of 1) and dense forest with scrub forest (resistance value of 1-20), including: percentage of habitat within the focal landscape (PLAND); patch density (PD), a measure of the subdivision of habitat within the landscape (i.e., number of patches divided by area); and largest patch index (LPI), the percentage of the focal landscape represented by the largest patch of a given habitat type parameters in our models were based on other studies on tigers and other wide-ranging felids (Thatte et al. 2018;Kaszta et al. 2019Kaszta et al. , 2020. Simulated movement across the landscape in this approach requires two input layers: (1) a cost-distance matrix; and (2) a set of initial population source locations. Cost-distance matrices, which define the cost-weighted distance between every pair of potentially occupiable locations, were generated in UNICOR (Landguth et al. 2012) for each combination of landscape change scenario (SCN), resistance surface transformation (RST), and population density (DEN). To produce initial source points for simulations, we selected a sample of 30 points, representing the upper range of the current tiger population estimates within DPKY boundaries (Ash et al. 2020b), generated probabilistically and proportional to predicted tiger habitat suitability (Ash et al. 2021). These source points served to represent the first generation in our simulations. As this is a non-overlapping generation model, following breeding, all adults are removed. Further mortality of individuals is determined by one of two processes: the absence of available territories within a cost-distance threshold and as a function of the probability of mortality (see "Mortality Functions"). Additional model details and a full summary of parameters can be found in Online Resource 1.

Mortality Functions
Following the simulation of individual movement and breeding within a single generation of our simulations, we applied a series of spatially-explicit mortality functions which define probability of mortality for each cell in our study extent via habitat suitability, human presence, and other conditions. The explicit incorporation of spatially-differential mortality risk in landscape-scale population modelling studies has been shown to dramatically affect conclusions pertaining to landscape connectivity and population persistence (Kramer-Schadt et al. 2004;Cushman 2006;Kaszta et al. 2019Kaszta et al. , 2020. We aimed to evaluate the effect of spatially-differential mortality risk (MRT) and sensitivity of simulations to different mortality function development approaches based on landscape resistance, predicted poacher presence, expert-based surfaces, and an effective population size model. We developed a total of 13 functions using these approaches, representing the probability of a tiger that has dispersed to a given pixel dying before reproducing. We present these in Table 3 and describe them in additional detail in Online Resource 2.
For comparative purposes, we tested a null mortality function in our simulations (M0) in which no additional mortality was applied. In effect, this function represents the maximum theoretical population growth rate, dictated by resistance surface, movement, available habitat, and demographic parameters.
We developed a habitat-based mortality approach in which mortality is a function of local, empiricallybased predicted tiger habitat suitability (MSSO model from Ash et al. 2021), landscape change scenario, resistance surface, and protected area coverage. This approach assumes mortality risk is higher where habitat-suitability is low and resistance is high, with greater risk outside protected areas. To test sensitivity, we transformed resulting mortality surfaces by a range of power functions (0.5, 0.75, 1.0, 1.5, and 2.0), rescaling each between 1-100 to reflect varying intensity and non-linear probabilities of mortality. A total of five functions were developed using this approach (MH05, MH075, MH10, MH15, MH20).
We tested an additional mortality function (MP1) derived from a local empirical model of predicted poaching presence. This was generated using a random forest model in which environmental and anthropogenic factors served as explanatory variables, with poacher presence generated from camera-trap surveys (Ash et al. 2020d) as the response, averaged across 10 bootstrapped iterations. In this function, predicted poaching presence (including both timber and wildlife poachers) is assumed to be correlated with overall probability of mortality (1-100) due to the risk of direct mortality, such as from targeted hunting/snaring, and indirect effects, such as from prey depletion.
We also incorporated four expert-derived functions (ME1, ME2, ME3, and ME4) based on functions tested in Ash et al. (2020a). This approach defines spatially-differential mortality risk based on available tree cover, protected area coverage, and road presence. In ME1 and ME2, the effect of roads was included via rescaled road density within a 20 km radius. In ME3 and ME4, the effect of roads was included as a function of distance.
Lastly, we also tested two additional mortality functions used in Ash et al. (2020a;MR1 and MR2), derived from an India-based effective tiger population model (Reddy et al. 2017). In this approach, probability of mortality is a function of broad-scale protected area coverage, mean forest cover, and mean landscape resistance with the assumption that mortality risk is inversely proportional to predicted effective population size.  (2005) Table 1).
Here, we discuss overall model sensitivity results with a detailed assessment of the effect of landscape change scenario and other factors under the same mortality function approach.

Model sensitivity
Simulated N 20 were most sensitive to mortality function (MRT), with dramatic differences in population trajectory based on mortality approach (Fig. 3,  Fig. 4; Online Resource 3, Fig. 1-6). Simulations with no additional mortality (M0) resulted in populations expanding to carrying capacity, dictated by density (DEN; 100% [n = 990] of simulations resulted in population growth; Table 4 Patterns of simulated N 20 when adjusting landscape change scenario (SCN) were not as strong as those when adjusting mortality, but were nonetheless clear ( Table 4). The percentage of simulations which resulted in extinction ranged from 61.5% (n = 720) for scenario R1 (habitat restoration) to 77% (n = 901) for scenario D0 (combined development). Conversely, the percentage of simulations which resulted in stable or increasing populations ranged from 14.7% (n = 172) for D0 to 30.4% (n = 356) for R0 (combined restoration). Predicted N 20 was less sensitive to nonlinear transformation of landscape resistance (RST). Extinction rates were generally lower in simulations with higher, concave landscape resistance. Specifically, simulations resulted in extinction in 64.5% (n = 2,765) of cases for RST05 and 70.9% (n = 3,040) of cases for the convex RST20. In comparison, the percentage of simulations which resulted in stable or increasing populations ranged from 22.4% (n = 960) under RST20 to 26.8% (n = 1,148) under RST05.

Results under moderate mortality
Mortality had a disproportionately large effect in our study as the majority of mortality functions produced extreme values (N 20 = 0 or N 20 = K), precluding an assessment of nuanced differences between factors such as landscape change scenario. Here, we review results holding mortality constant at moderate levels (MH10), under which N 20 exhibited particular sensitivity to parametric variation.

Landscape change scenario
Holding mortality function constant (MH10), N 20 was most sensitive to differences in landscape change scenario (Table 5; Online Resource 3, Table 2). Under the base scenario (B0), median population at generation 20 ( Ñ 20 ) was 27.5, differing only slightly from starting N 1 (30). Population trajectories were varied-approximately 48.9% of simulations resulted in stable or increasing populations, 27.8% resulted in population declines, and 23.3% ended in extinction.
Development had an almost universally negative influence on population trajectory. A substantial majority of simulations with development-based landscape change scenarios resulted in population declines or extinctions. When comparing results between landscape change scenario and B0, half of scenarios resulted in very highly significant differences in mean population ( N 20 ; p < 0.0001) while half resulted in no significant difference (Table 5).
The combined development and restoration scenario (DR0), produced the greatest difference in observed mean of N 20 compared to the base scenario ( N 20 = − 39.9***) and had an Ñ 20 of 0. Approximately 95.6% of simulations under this scenario resulted in extinction, 4.4% resulted in population declines, and 0% produced population increases. Results under the combined development scenario (D0) were similarly severe. Under D0, simulations had a mean difference in N 20 of − 38.3 (***) compared to B0, Ñ 20 = 0, and extinctions or declines in 83.3% and 15.6% of cases, respectively. The landscape change scenario where overall road use/resistance increased (D2) was also notably severe. Scenario D2 had a mean difference of N 20 = − 35.8(***) compared to B0, Ñ 20 = 0, and extinction and population declines in 74.4% and 20.0% of cases, respectively. The landscape change scenario incorporating the construction of additional dams/reservoirs (D4) also had a comparatively negative effect ( N 20 = − 27.0***, Ñ 20 = 0); in D4, 52.2% of simulations resulted in extinction, 30.0% resulted in population declines, and 17.8% resulted in stable/increasing populations Development of Route 3462 through DPKY (D3) also produced significantly lower mean population ( N 20 = − 22.9***) compared to B0, and low relative Ñ 20 (9.5). Under this scenario, approximately 45.6% of simulations resulted in extinction, 32.2% resulted in population declines, and 22.2% resulted in stable or increasing populations. These development scenarios were similarly strong in their effect, with differences in N 20 only significant between D3 and DR0/D0(*).
Significant differences in N 20 were not observed for two development scenarios: D1, in which tree cover outside PAs was converted to village/agriculture matrix; and D5, representing doubling of urban area coverage and development of border crossings. Notably, in all restoration scenarios tested (R1 -reforestation/increased PA coverage; R2 -wildlife crossing structures; R0 -all restoration measures), N 20 did not differ significantly from the base scenario (B0).

Resistance transformation
Varying resistance surface transformation (RST), holding mortality constant at MH10, had similar, clear patterns of difference in N 20 as observed in overall results (Tables 4 and 5). Concave/higher resistance (RST05) produced an N 20 = 28.8, and simulations resulted in extinctions in 37.6% of cases, declines in 26.4% of cases, and stable/increasing populations in 36.1% of cases. In contrast, convex/ lower resistance (RST20) had an N 20 = 21.2 and simulations resulted in extinction in 53.0% of cases, percentage of population declines in 22.4% of cases, and stable/increasing populations in 24.5% of cases. Differences in observed N 20 were only significant between RST20-RST05 ( − 7.6*; Online Resource 3, Table 4).

Interactions between factors
Results from ANOVA under moderate mortality (MH10) were significant for all main effects (Online Resource 3, Table 2), while two-way and three-way interactions between factors were not. When examining specific interactions within and between factor groups with Tukey HSD testing, the majority of significant differences were observed when landscape change scenario (SCN) varied. Only one significant interaction was observed when holding SCN constant (B0:6KM:RST20-B0:6KM:RST05; N 20 = -61.2*; Online Resource 3, Table 5). There were only five instances where interactions with other factors produced significant differences in N 20 between B0 and landscape change scenarios not considered significant without interactions (

Discussion
Population dynamics are driven by demographic processes interacting with the influences of heterogeneous landscapes. In this study, we simulated population dynamics of an isolated tiger population of conservation priority using an individual-based, spatially-explicit approach. Our goal was to quantify the relative effects of landscape change scenarios on population persistence as well as sensitivity of predictions to spatially-explicit mortality risk, resistance to dispersal, and population density. Our most important finding is that spatially-differential mortality risk, and the approach used to define this risk, dominated predictions of population persistence. Second, we documented significant differences in population trajectories depending on landscape change scenario. Notably, we documented strong negative effects of development on population size and persistence, even when habitat restoration measures were applied. Resistance transformation and its interaction with mortality had a marginal effect, while adjusting maximum population density did not substantially influence population trajectory. The high sensitivity of predictions to spatially-explicit patterns of mortality risk, and the high variability in simulated population trajectories among model runs, starkly highlight the vulnerability of this small and isolated, but regionally important tiger population. Our results provide clear predictions to aid in the development of management and conservation strategies for this population, as well as spatiallyexplicit population models of other threatened species. In addition, our study also provides valuable insight to inform the development of spatially-explicit population models for other threatened species. Our results demonstrate that management efforts to increase connectivity is not likely to be beneficial if it is not coupled with efforts to mitigate mortality risk, and indeed, increased connectivity can increase extinction risk in situations where mortality risk is not mitigated (e.g., Cushman 2006;Cushman et al. 2010).

Mortality
Predictions of population trajectory and persistence were overwhelmingly driven by spatially-differential mortality risk. Among the 13 mortality functions applied across our simulations, nine produced extinction rates greater than 50% with five resulting in extinction rates of 100% by generation 20. Concurrently, we documented significant differences in predicted population depending on the mortality approach applied, with predictions notably sensitive to minor adjustments in probability of mortality. Moderate, curvilinear adjustments in habitatbased mortality, for example, produced opposing predictions of population trajectory. Overall, this strong sensitivity to mortality risk and the observation of frequent or universal extinction among realistic mortality functions highlights, (1) the high risk faced by this small and isolated tiger population, and (2) the dominant effect of mortality risk in driving population dynamics and extinction risk of small populations.
The clear and strong effect of mortality in our simulations has important implications. Our results are consistent with a related population simulation study conducted on tigers in this landscape in which the incorporation of spatially-differential mortality risk dramatically shifted predictions of population dynamics, extinction risk, and broad-scale landscape connectivity (Ash et al. 2020a). Specifically, Ash et al. (2020a) used a different modelling approach, employing spatially-dynamic resistance kernel modelling (Cushman 2015;Barros et al. 2019) to evaluate population size, extinction risk and spread. Their study found that, under all realistic scenarios of mortality risk, population range expansion was highly constrained and population declines or extinctions were frequent. Other tiger population modelling studies have reported that marginal increases in poaching rates produced large increases in predicted extinction risk (Kenney et al. 1995;Linkie et al. 2006;Chapron et al. 2008) and reductions in population connectivity (Carroll and Miquelle 2006), particularly for small, isolated populations. These trends have also been observed in simulations of other large felids (Kramer-Schadt et al. 2004;Newby et al. 2013;Kaszta et al. 2019;Naude et al. 2020). It is evident from our simulations that the degree of mortality risk within DPKY is likely the most critical determining factor for the future of this population. Thus, detecting and mitigating sources of mortality will likely be critical to achieving short-and long-term tiger population persistence in this landscape. Importantly, efforts to increase population density (such as by improving habitat or prey density) or improve movement across the landscape could be destined to fail if mortality risk is not effectively addressed.
The disproportionate effect of mortality in our study also has implications for conservation of other small populations and the development of similar spatially-explicit population models for other species. Broadly, our results appear to reflect the inherent vulnerability of small populations to stochastic processes, allee effects, and related dynamics of extinction risk (Caughley 1994;Dennis 2002;Fagan and Holmes 2006). In such small populations, even marginal increases in mortality risk may act in concert with demographic stochasticity to drive populations to extinction. More specifically, our findings are similar to other modelling studies in which the inclusion of spatially-differential mortality risk dramatically changed predictions of population dynamics, connectivity, and the effect of large-scale development (e.g., Kramer-Schadt et al. 2004;Cushman 2006;Kaszta et al. 2019Kaszta et al. , 2020. Kaszta et al. (2019) and Kramer-Schadt et al. (2004) highlight that the deleterious impacts of development/infrastructure are evident not only in their effect on habitat, but also in reducing population viability through reduced connectivity and road deaths. Further, as evidenced in our study and others (Kaszta et al. 2019;Ash et al. 2020a), differences in the approach used to parameterize mortality across space and its degree of severity, even to a marginal degree, can have considerable influence on model predictions. This underscores the importance not only of incorporating spatial variation in mortality risk in population simulations, but conducting sensitivity analysis to evaluate its degree of influence on predictions. For similar studies, the lack of inclusion of spatial mortality risk, or parameterization that neglects key factors, can result in the underestimation of the effects of infrastructure and undermine prediction accuracy.

Landscape change scenarios
A central goal of our study was to discern the relative effects of various landscape change scenarios on predicted long-term population dynamics. Notably, we documented strong and significant differences in predicted populations between current conditions and scenarios of landscape change, when holding mortality risk constant at moderate levels.
We documented strong negative effects in scenarios which simulated potential expansion of infrastructure development. Simulations in which the primary change to the landscape involved intensification or development of major roadways in the DPKY forest complex (D2/D3) resulted in significantly lower predicted population size, with relatively high predicted rates of extinction. Roads are commonly associated with negative effects on biodiversity (Coffin 2007). These effects have been known to exert a particularly strong influence on the viability and connectivity of populations of tigers (Kerley et al. 2002;Thatte et al. 2018) and other large carnivores (Kramer-Schadt et al. 2004;Kaszta et al. 2019). In recent studies on clouded leopards (Neofelis nebulosa and N. diardi), which applied a similar spatially-explicit approach to evaluate relative effects of development scenarios, linear structures such as roads and railways were associated with significant declines in predicted population size (Kaszta et al. 2019(Kaszta et al. , 2020. Notably, our results are consistent with other evidence from the DPKY complex which documented broad-scale, negative associations between tigers and roads (Ash et al. 2020d).
Development of dams/reservoirs along the periphery of protected areas had a similarly strong and significant negative effect. While a number of studies have been dedicated to quantifying the varied ecological impacts of dams (Jones et al. 2016;Wu et al. 2019), investigations into the potential implications for terrestrial mammals upstream from such projects are comparatively lacking. The impact of dams for wide-ranging carnivores may be greater than localized habitat loss (Kaszta et al. 2020). It is possible that reservoirs created by dams may improve access into remote areas for illicit activities such as poaching in a similar way to roads (Bennett and Robinson 2000) and could catalyse further infrastructure development or land use change in their vicinity (Lillesund et al. 2017). For small populations such as DPKY, even marginal associated increases in mortality risk and disturbance along park peripheries could be enough to shift population trajectories.
We did not find significant differences between predicted population in the base scenario and scenarios simulating habitat loss outside protected areas (D1) or expansion of urban/border areas (D5). It is possible that this is due to the primary area of change occurring within an already human-dominated matrix of high resistance and mortality risk. Importantly, the structure of the DPKY system is represented by a single block of mid-to high-quality tiger habitat (Ash et al. 2021) that is not substantially subdivided by a developed matrix. Thus, the fact that tigers do not reside in, nor are required to move through the surrounding matrix, may explain the lack of strong effects of additional development outside this block. For areas in which tigers occur in a metapopulation, conditions in the matrix between populations could have considerable implications for movement and persistence (Thatte et al. 2018).
Restoration scenarios such as reforestation with expanded protected area coverage (R1) and installation of highway crossing structures did not significantly differ in associated predictions compared to the base scenario. This may be attributed to the smallscale nature of these simulated changes relative to the size of the landscape. In the case of highway crossing structures, their potential effect on simulated population from connectivity may be marginal relative to mortality risk, which was not measurably reduced by these structures and was the dominant driver of population dynamics in our study. Our results are consistent with a related study in DPKY which did not document significant effects of such structures on simulated population and connectivity, particularly compared with other factors such as mortality (Ash et al. 2020a).
Importantly, the impacts of landscape change scenarios in our simulation extend beyond physical changes to the landscape, but also include the effect of these changes on the spatial pattern of mortality risk. Predicted population dynamics in our study is strongly influenced by synergistic effects between landscape patterns and mortality risk. Studies of threatened populations which do not account for the associated variation in mortality risk with landscape change may underestimate the effect of these changes on population persistence (Thatte et al. 2018;Kaszta et al. 2019). Furthermore, the strong effect of landscape change in our study and inseparable link between demography and landscape dynamics (Lande 1988;Kareiva and Wennergren 1995) reinforces the importance of considering future potential landscape configurations in long-term population modelling. This is particularly important in terms of how landscape change affects spatial patterns of mortality risk.

Combined landscape change scenarios
In addition to the independent evaluation of specific landscape changes, we also explored differences in predicted population when these changes were combined. Scenarios which incorporated all development measures (D0 and DR0) had a profoundly negative effect on predicted population size and produced very high extinction rates. Importantly, incorporation of restoration measures failed to even marginally mitigate the simulated negative effects of extensive development and its associated increases in mortality risk.
Where infrastructure is likely to negatively impact tiger populations, measures such as road crossing structures for wildlife, rehabilitation and restoration of habitat, or financial compensation have been recommended as a means to offset their impacts (Quintero et al. 2010). However, a degree of caution is merited when considering such offsets. Results from our study and others suggest that these measures may be insufficient in ameliorating the long-term implications of habitat loss, fragmentation, and increased mortality risk resulting from development scenarios and may fail to justify new infrastructure development in the habitat of vulnerable species (Ash et al. 2020a). Rehabilitation and restoration of habitat must result from strategic planning and evaluation, and must be designed for specific purposes such as enhancing or enabling population connectivity. For example, simulations by Kaszta et al. (2019) suggested that targeted restoration of forest corridors could mitigate the effects of road development on clouded leopard population size and connectivity to a limited degree. However, the same study documented an interactive effect between the re-alignment of certain linear structures and increased mortality, depressing population size relative to original plans. Both our study and Kaszta et al. (2019) underscore the potentially catastrophic synergy between development and spatial patterns of mortality risk in driving patterns of extinction.

Resistance transformation
While sensitivity analysis suggested that landscape resistance to dispersal had lower comparative influence on model predictions, we documented a significant difference between predicted population size based on the form and intensity of landscape resistance. Higher, concave resistance (RST05) resulted in higher predicted mean population size relative to lower, convex resistance (RST20) and had an interactive effect with mortality. High resistance may preclude movement critical to persistence in a metapopulation (Cushman 2006;Rudnick et al. 2012). However, movement by tigers into human-dominated matrices is often highly risky due to the risk of direct mortality, conflict with humans, and lack of natural prey (Smith 1993;Goodrich et al. 2008). In our isolated population, high resistance may have had the effect of discouraging and reducing movement to areas of high-risk that would otherwise increase mortality probability. Specifically, when landscape resistance was high outside of optimal, low-risk habitat, populations were predicted to be greater because dispersing individuals avoided high-risk areas due to increased resistance to movement. The ability of a vagile species to navigate an unfriendly landscape may increase the effect of landscape barriers and dispersal-related mortality on population viability relative to those with limited movement ability (Carr and Fahrig 2001;Carr et al. 2002;Kramer-Schadt et al. 2004;Ash et al. 2020a). It is possible, albeit perhaps non-intuitive, that a highly disturbed matrix could yield benefits to certain species compared to one that is partially disrupted, if it discourages movement through or settlement in areas with low suitability and high risk (Cantrell and Cosner 1993). This may also be the case if such barriers inhibit movement to ecological traps, such attractive sink habitat, that may provide potential suitability, but with high mortality risk (Woodroffe and Ginsberg 1998;Delibes et al. 2001). Further, if negative effects of movement in such areas exceed potential benefits, this may ultimately undermine population persistence (Burkey 1989). Such patterns have been observed in studies of other species which have reported high dispersal associated with declining populations where mortality risk is high (Cushman 2006;Cushman et al. 2010). However, studies on the real-world effects of these patterns on large carnivore movement and persistence are lacking (Cushman et al. 2016).

Population density
Adjustment of maximum population density in our study had a minimal effect, primarily dictating the ceiling of simulated population in simulations incorporating little to no additional mortality risk. A study by Karanth and Stith (1999) suggests that tiger populations with a higher carrying capacity, supported by high prey abundance, may be able to persist despite elevated mortality, such as from poaching. We did not detect a clear pattern in which high population carrying capacity measurably reduced negative impacts on simulated population from mortality or landscape change scenario. This reinforces our findings of mortality acting as the primary driver of population dynamics. Further, these results are consistent with other studies highlighting the potentially severe, longterm effect of increases in mortality on population persistence, even in relatively large populations or those with high potential density (Kenney et al. 1995;Chapron et al. 2008;Tian et al. 2011).

Limitations/caveats
The strength of our assessment lies in the comparison of the relative effects and interactions of factors across a realistic range of landscape change scenarios and parameters in order to inform potential management interventions and areas of future research (Dunning Jr et al. 1995;Reed et al. 2002;Linkie et al. 2006;Kaszta et al. 2019). Our approach of incorporating sensitivity analysis is particularly useful given the inherent uncertainty in defining parameters for such simulations (Boyce 1992;Dunning Jr et al. 1995;Wiegand et al. 2004). Where available, our individual-based, spatially-explicit simulations are guided by empirical data (Ash et al. 2020b(Ash et al. , c, 2021 and credible parameters governing tiger behaviour utilized in similar simulations elsewhere in the tiger's range (Thatte et al. 2018). This robust approach grounds our assessment within the realm of plausible reality. However, we wish to underscore that the focus of this study was not to make predictions of specific population sizes in relation to landscape change scenarios or mortality rates. Our findings are beneficial in identifying the types of landscape change likely to be most impactful among those tested, and degree of sensitivity of this population to key factors (notably mortality). These insights may aid in official planning where they can be compared with specific development and management plans in concert with up-todate information on this tiger population and factors that could affect population persistence.
We note that exploration of restorative landscape changes or assessment of alternative development scenarios (such as re-alignment of infrastructure; sensu Kaszta et al. [2019]) are relatively limited. Given the high degree of long-established human presence outside protected area boundaries, opportunities to evaluate broader-scale habitat restoration and other measures were inherently constrained in this landscape. We did not investigate long-term genetic trajectories in this population, though this too will likely be a factor in population persistence, which can be severely undermined by inbreeding depression (van Noordwijk 1994;Vasudev et al. 2017). Notably, while mortality probability varied across space, it was held constant over our simulations; in reality, there is certain to be stochastic variation in this and other factors over time. We also did not consider potential catastrophic environmental stochasticity, such as disease, whose effects may be further exacerbated in small populations (Harihar et al. 2018). The potentially strong and pervasive effects of climate change may also threaten population persistence across landscapes (Marks 2011;Rudnick et al. 2012;Wasserman et al. 2012;Dar et al. 2021). We recommend additional, targeted studies to evaluate the impacts of potential associated habitat changes. However, it remains clear from our study that such long-term factors for tigers in this landscape may be moot if the calamitous effects of elevated mortality risk are ignored.

Conservation implications/recommendations
The tiger population in DPKY has been described as a population of national, regional, and global priority (Ash et al. 2020c), particularly given the catastrophic range and population declines elsewhere in Southeast Asia (Sanderson et al. 2006;Lynam and Nowell 2011). Our results suggest that this population may be situated on a proverbial knife's edge. Predictions of population trajectory and persistence in our study were highly sensitive, particularly to mortality risk. Importantly, we found that marginal shifts in mortality, particularly combined with landscape change, drastically shifted predictions. When holding mortality risk constant at moderate levels (MH10), roughly half of simulations produced extinctions. When transforming this function to slightly convex (MH15) and concave (MH075) forms, extinction rates diverged to 12% to 85%, respectively. These results reflect the high degree of sensitivity and vulnerability of small, isolated populations to stochastic processes (Kenney et al. 1995;Vasudev et al. 2017;Thatte et al. 2018).
The tiger population in DPKY is not unique in respect to being smaller than recommended for longterm viability (Wikramanayake et al. 2011;Kenney et al. 2014). Such populations may be supported by rescue effects through immigration (Linkie et al. 2006;Banerjee et al. 2010;Thatte et al. 2018). However, many populations, including DPKY, are likely too isolated for this to be the case without active intervention. Small populations may benefit from the tiger's naturally high fecundity (Karanth and Stith 1999), though high mortality, combined with isolation, can quickly drive populations to extinction (Chapron et al. 2008). Empty forests throughout Southeast Asia (Lynam and Nowell 2011;Rasphone et al. 2019;Ash et al. 2020c) highlight that requisite factors such as presence of forest cover and potential habitat mean little for species like tigers if unsustainable mortality is not prevented.
Evidence from our study, consistent with others on tigers (Kenney et al. 1995;Linkie et al. 2006;Chapron et al. 2008), reinforces the point that preventing mortality, whether from poaching or development in tiger habitat, is likely the most critical measure to be considered in developing management and conservation strategies.
Infrastructure development, specifically expansion and intensification of roads and dam construction, were associated with profoundly negative effects on predicted tiger population persistence and was notably synergistic with elevated mortality. Our results strongly reinforce recommendations in Thailand's national tiger action plan against the development of infrastructure in tiger habitat (Pisdamkam et al. 2010). Expansion of highways such as Route 3462 and construction of dams within protected areas (e.g., Huay Satone) are incompatible with this recommendation. Predicted negative effects of dams also merits targeted research on the potential effects of dams on tigers and other terrestrial mammals. Further, efforts to mitigate the effects of infrastructure, such as wildlife crossing structures along roads should be planned carefully, given that use of these structures may vary due to complex factors and by species (Clevenger and Waltho 2000;Caldwell and Klip 2020). Specific research on the use and effect of current structures on tigers is also recommended.
Given the vulnerability of this population, our results reinforce previous recommendations (Ash et al. 2021) for the prioritization of core habitat in central DPKY as an area of critical conservation priority. Efforts to conserve and effectively manage this population should be underpinned by highstandard and regular monitoring of the tiger population, prey, and threats. Designing comprehensive and timely population monitoring and enforcement patrolling is critical in order to detect and respond to emerging threats such as poaching (Duangchantrasiri et al. 2016) or disease (e.g., canine distemper virus, Seimon et al. 2013;Terio and Craft 2013), which could quickly overwhelm small, vulnerable populations. DPKY's recently established Dong Phayayen-Khao Yai Wildlife Research Station could serve as an effective foundation for these efforts, especially for providing a centralized, landscapewide monitoring location and acting as a local training centre for methods of wildlife research.
For isolated populations such as DPKY, active population management measures, such as translocations, may be an option to prevent the effects from inbreeding depression from being realized over longer time periods (Kenney et al. 2014). Additional simulations investigating the potential effect of translocations on heterozygosity and allelic richness may be particularly helpful in informing the development of such measures. Lastly, we believe our approach can be helpful for investigating the effects of dynamic and changing landscapes on other populations of tigers at similar conservation crossroads.