Occupancy‐based monitoring of ungulate prey species in Thailand indicates population stability, but limited recovery

Longitudinal studies of wildlife are urgently needed in South-East Asia to understand population responses to the high poaching pressure that characterizes this region. We monitored population trends and habitat use of five heavily poached ungulate species (gaur, sambar, wild pig, red muntjac, and Fea’s muntjac) over five years in two protected areas in western Thailand using camera trap surveys. We used single-season occupancy models to investigate effects of ecological and anthropological variables on ungulate distribution, and multi-season models to assess occupancy dynamics over time. Occupancy of gaur and sambar was low (<0.25), but concentrated near saltlicks and at low elevations. Wild pig and muntjac occupancies were 3–4 times higher (0.60–0.80). Wild pig occupancy was lower near villages, but this effect dissipated in the final year of the study, coinciding with a purported decrease in poaching. Wild pig occupancy increased significantly, with the probability of colonizing new sites doubling from 0.40 to 0.81 over time. In contrast, occupancy rates of gaur, sambar, and muntjac did not grow, though they were stable. Poaching pressure during the study was low, perhaps allowing populations to stabilize. But only wild pig (the most resilient of the five species) increased. The failure of gaur and sambar to recover might stem from historical overhunting combined with ecological constraints, particularly low saltlick density. Recovery of ungulates (and the carnivores that depend on them) in overhunted South-East Asian reserves might require intensive interventions, particularly habitat improvement and population augmentation, to achieve conservation objectives.


INTRODUCTION
The conservation status of ungulates and carnivores is declining more rapidly in South-East Asia than any other region of the world due to habitat loss and overhunting (Di Marco et al. 2014). Population declines and rising extinction risks are particularly severe for large-bodied carnivore and herbivore species over 100 kg (Ripple et al. 2014(Ripple et al. , 2015. Ungulates such as deer and wild cattle are key prey items for endangered large carnivores of Asia (Hayward et al. 2012), influencing carnivore habitat use, density, reproduction, competitive interactions, and prospects for coexistence Stith 1999, Steinmetz et al. 2013). The depletion of ungulate prey is thus a major threat to the survival of the carnivores that depend on them (Wolf and Ripple 2016). For example, prey depletion reduces reproduction and lowers habitat carrying capacity for tigers (Panthera tigris; Karanth and Stith 1999), hindering tiger population recovery (Sanderson et al. 2006). A major conservation challenge, globally and especially in South-East Asia, is to protect and recover ungulate communities, for their own survival as well as the carnivore communities that rely on them.
Protected areas in South-East Asia often have sparse ungulate populations that have been reduced from illegal hunting (i.e., poaching). Management interventions such as ranger patrolling are often used to reduce pressure on wildlife and thereby instigate their recovery. A key task of protected area management is monitoring the status of its wildlife populations over time in response to management. The influence of ecological aspects must also be considered when assessing the status of ungulates. For example, habitat preferences and food availability could affect patterns of population recovery independent of human pressures (Steinmetz et al. 2010). Longitudinal studies (i.e., over time) of wildlife are uncommon in South-East Asia compared to other regions, yet are urgently needed to understand population dynamics with respect to the intense poaching pressure faced by many species here (de Silva 2016).
In this study, we monitored the population status of five ungulate species recovering from overhunting in two protected areas in Thailand. We investigated the effects of ecological and anthropogenic variables on the distribution and population trends of gaur (Bos gaurus), sambar (Rusa unicolor), wild pig (Sus scrofa), northern red muntjac (Muntiacus vaginalis), and Fea's muntjac (M. feae). These ungulates constitute key prey species for three large carnivore species that inhabit the region: tiger, leopard (P. pardus), and dhole (Cuon alpinus; Sunquist 1995, Simcharoen et al. 2018). Gaur (650-900 kg) and sambar (185-260 kg) are vulnerable to global extinction, with commercial poaching driving declining populations (Duckworth et al. 2016, Timmins et al. 2016a. Red muntjac (20-28 kg) is still widespread in South-East Asia though many populations are declining (Timmins et al. 2016c). Fea's muntjac is limited to a small portion of Thailand and Myanmar, and little is known of its status (Timmins et al. 2016b). Wild pigs (75-200 kg) are widespread and abundant in many places, though some populations are reduced where poaching is intensive (Oliver and Leus 2008).
Ungulates were difficult to observe directly due to their low densities, often nocturnal habits, and low visibility inherent to their dense forest habitats. These conditions hindered the use of methods to estimate animal abundance directly, such as distance sampling. Instead, we used occupancy as the state variable to monitor population change. Occupancy models account for imperfect detection of the target species, yielding estimates of probability of occurrence and influence of variables on occurrence. For rare species that are difficult to observe, occupancy-based monitoring programs offer a valuable alternative to methods that aim to track population size directly (Guillera-Arroita et al. 2010) and have been employed to monitor reptiles (Sewell et al. 2010), birds (Olson et al. 2005), and mammals (Steinmetz et al. 2014). Our aims were to understand (1) the drivers of ungulate distribution in the landscape, (2) whether ungulate populations were growing, and (3) if so, which variables affected population growth. For the first question, we used single-season occupancy modeling. For the latter two questions, we used multi-season occupancy modeling.

Study sites
The study was conducted in Mae Wong National Park (894 km 2 ) and Khlong Lan National Park (300 km 2 ) in the Dawna Mountain Range, western Thailand (99°4 0 30″ E to 99°22 0 13″ E and 15°39 0 52″ N to 16°20 0 35″ N; Fig. 1). Elevations range from 150 to 1964 m above sea level. The area has a monsoonal climate with a dry season (November-May) and a wet season (May-October). Average annual rainfall is 1200 mm. Mean temperature is 27°C. The major habitat types are mixed deciduous forest, secondary forest with bamboo, and semi-evergreen forest.
Both parks were subject historically to intensive commercial resource use. Logging occurred throughout most of the parks from the 1970s v www.esajournals.org until 1985 in Khlong Lan and 1987 in Mae Wong (when the parks were first established). This was accompanied by commercial wildlife poaching to supply markets in surrounding towns (WEF-COM 2004). Populations of ungulates plummeted during this era (Wildlife Conservation Bureau 1994, Pattanavibool andPhoonjampa 2017), and the parks have depressed densities of ungulates and large carnivores.
Mae Wong and Khlong Lan national parks (hereafter, MWKL) are part of Thailand's Western Forest Complex, a complex of protected areas that covers 19,000 km 2 and harbors one of the largest populations of tigers (Simcharoen et al. 2007) and probably leopards and dholes (Kamler et al. 2015, Rostro-García et al. 2016 remaining in mainland South-East Asia. Tiger density is a function of prey density and biomass (Karanth et al. 2004). Less than 4000 wild tigers remain in Asia, and the 13 tiger-range countries have committed to doubling the global tiger population (Global Tiger Initiative Secretariat. 2011). Sites like MWKL, where tiger populations are well below potential carrying capacity, could have a major role in this recovery, but only if prey populations can be increased (Harihar et al. 2018). Thus, the recovery of ungulate prey species here, and the lessons learned while attempting to do so, is of global importance for recovery of tigers.

Camera trap surveys
Passive infrared camera traps (Cuddeback, Attack IR, Ambush, and C1 models, Cuddeback Company; and Bushnell, Trophy Cam and Trophy Cam Aggressor models, Bushnell Corporation) were set across the two parks in 2012, 2014, and 2016 as part of a long-term tiger and prey monitoring effort (Fig. 2). We set cameras in pairs, with one camera set to take photographs, and the opposing one taking a 20 s video. We established 67 camera sites in the first year. In v www.esajournals.org subsequent years, most camera sites (80-82%) were at the same location, or within 1500 m of the original location. A minority of sites (18-20%) were different between years, but the same general areas were sampled every year (Fig. 2), and with similar effort (4.9-5.2 camera locations/100 km 2 each year); thus, the same ungulate populations would have been sampled over time despite some camera sites shifting location. We sampled each site for 8-12 weeks and used 1-week periods as occasions for occupancy modeling. We did not have sufficient staff or cameras to cover the entire study area at once, so we divided the area into three equal-sized zones and sequentially sampled each zone for the desired period. The initial zone was rotated each year. Sampling was mainly in the dry season each year, starting in December (access is reduced in the rainy season due to flooding), so inferences about habitat use pertain specifically to the dry season.
Camera traps were attached to trees at approximately 50 cm above the ground and set to operate day and night, with a 30-s interval between captures. Average spacing between cameras was 2.6 km. Cameras were placed along ridges, streams, animal trails, and old roads, which are landscape features used by the focal species. Separating the two muntjac species from black and white nighttime photos was difficult, so we lumped these species for analysis. Whereas muntjac and sambar are often solitary, gaur and wild pig populations are comprised of both solitary animals (typically single males) and herds. We lumped gaur and wild pig social units (single males, herds) together for occupancy analysis, as our goal was to investigate species-level habitat use and population trends. This approach incurs the possibility that herd size might increase through time while occupancy remained unchanged, obscuring our ability to discern important population changes. We examined herd sizes to account for this possibility. We counted numbers of animals from the 20-s videos and used these data to estimate herd size in each year. Animals usually walked in v www.esajournals.org one direction (not back and forth), which helped avoid double-counting. These counts are probably underestimates, but this bias would be consistent throughout the study. We examined differences in mean herd size across the three sampling periods using Kruskal-Wallis tests.

Ecological and anthropogenic variables
We modeled occupancy as a function of five variables that could plausibly influence ungulate distribution and recovery: (1) distance to nearest permanent stream, (2) distance to nearest natural saltlick, (3) distance to nearest village, (4) elevation, and (5) Normalized Difference Vegetation Index (NDVI) (further details in Appendix S1: Table S1). We expected that all species might occur closer to streams (our sampling was in the dry season), closer to natural saltlicks (which all five species use), and further from villages (poaching pressure is higher closer to villages in MWKL). We expected gaur and sambar, which are mixed grazers and browsers, to have higher occupancy at lower elevations where grass is more plentiful (Steinmetz et al. 2008, Simcharoen et al. 2014. NDVI is a satellite-derived surrogate for vegetation productivity, nitrogen content, and other aspects of high-quality food for herbivores, and is widely used to predict herbivore distribution and abundance (Pettorelli et al. 2011). We did not have an expectation for the effect of NDVI on our focal species, as they are habitat generalists (Bhattarai and Kindlmann 2012), but we included it as a potential explanatory variable that might help explain ungulate distributions. None of the variables was highly correlated with another (r < 0.56, Appendix S1: Table S2).
Distance to nearest stream, distance to nearest village, and elevation were obtained using a GIS database of the study area. Saltlick locations were recorded during field surveys and from the knowledge of park rangers. We used the MODIS (MODerate-resolution Imaging Spectroradiometer) MOD13Q1 V.5 dataset gathered at 250 × 250 m resolution, choosing data that corresponded spatially (pixels that overlapped with camera sites) and temporally (months and years that overlapped with camera operation) to our sampling. NDVI data were downloaded from the U.S. Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center (http://www.earthexploere.usgs.gov).

Single-season occupancy analysis
We used single-season, single-species occupancy models (MacKenzie et al. 2018) to identify ecological and anthropological factors underlying ungulate distributions. We built models for each species in each sampling year. Modeling was done in two steps. We first held occupancy constant while comparing two detection models: detection constant across sampling occasions (p(.)) and detection influenced by number of occasions that a camera was active (p(effort)). Retaining the best supported detection model as the base model, we then compared single-covariate occupancy models using AIC. If multiple variables were included in this top set, we created and assessed additive models as well. We looked for patterns of support for the same variables across years. We examined the beta coefficients (β) of variables in top models each year and their 90% confidence intervals; CIs that did not cross zero were considered statistically significant. We used liberal CIs in single-species occupancy analyses to avoid missing potentially important factors driving ungulate distributions; we expected traditional 95% CIs to be excessively wide for evaluating habitat covariate effects given the small sample sizes we had (Forthofer et al. 2007).
Occupancy models assume that sites (camera trap locations, in our case) are closed to changes in the state of occupancy during sampling. Our sampling occurred over 2-3 months, a short time relative to the life spans of the focal species, so this assumption was likely met with regard to demography (births, deaths). Camera spacing in our study (mean 2.6 km) was greater than the radius of home ranges of muntjak and sambar (Sukmasuang 2001, Leslie 2011) so for these species geographic closure was likely met as well. Dry season home ranges of gaur and wild pig are larger (Prayurasiddhi 1997) and might have overlapped with >1 camera, violating closure, so the occupancy estimator for those species should be interpreted as site use rather than occupancy (MacKenzie et al. 2002).

Multi-season occupancy analysis
Multi-season occupancy models estimate probabilities of occupancy, local colonization, and extinction. The colonization parameter, γ i (gamma), is the probability that a site previously unoccupied by the species in season i becomes v www.esajournals.org occupied in the subsequent season i + 1. The extinction parameter, ϵ i (epsilon), is the probability of an occupied site becoming unoccupied between seasons i and i + 1 (MacKenzie et al. 2018). Multi-season occupancy models assume populations are closed to changes in occupancy within seasons but allow changes between seasons. They are thus useful for monitoring changes in occupancy status over time and, additionally, can identify factors associated with trends in occupancy. The closure assumption applies within each season, similar to single-season models, but additionally it is assumed that there is no unmodeled heterogeneity in the colonization and extinction parameters. We met this assumption by modeling these parameters as functions of biologically plausible covariates.
We employed multi-season occupancy modeling to (1) estimate seasonal occupancy rates, (2) assess if occupancy increased over time, and, if so, (3) model the effect of covariates on that trend. We assessed detection probability first, comparing three detection models for each species (while holding occupancy and colonization constant): (1) detection constant across sample replicates and across years (p(.)), (2) detection varying by camera survey effort each year (p(effort)), and (3) detection varying across years (p(year)). We retained the best detection model for subsequent analyses.
Because we were interested in occupancy changes over time, we used a model parameterization which estimated seasonal occupancy, colonization rate, and detection probability (MacKenzie et al. 2018). We modeled occupancy as a linear function of time, treating occupancy in year 1 as a baseline, while allowing occupancy rates in subsequent years to follow a trend: logit(ψ) = β0 + β1(time). The design matrix for this model is in Appendix S1 (Table S3). We refer to this as the ψ(trend) model. The beta coefficient (β1) estimated by this model directly assesses population trend (on the logit scale). If β1 and its 95% confidence intervals were >0, we concluded that occupancy had increased for that species. This model also provided estimates of occupancy in each year, which we plotted to illustrate occupancy trends over time. We also estimated the rate of change in occupancy between successive years (λ t = ψ t+1 /ψ t ), which is analogous to growth rate of a population (MacKenzie et al. 2018). We looked for rates that were significantly larger than 1, indicating an increase in occupancy. For multi-season analyses, we used 95% confidence intervals to assess effect sizes (in contrast to 90% CIs in single-season analyses), in order to be highly confident in inferring that trends in occupancy rates of ungulate populations had occurred.
The colonization parameter was initially held constant as the focus was on change in occupancy, not underlying processes. If we detected a trend, then we modeled the effect of covariates on colonization, using the same covariates as in single-species modeling. We used program PRE-SENCE version 12.0 (Hines 2006) for all analyses. Models were ranked using Akaike's Information Criterion (AIC), and we selected the best model based on lowest AIC scores. We considered all models with delta AIC < 2 to have substantial support (Burnham and Anderson 2002).
By choosing occupancy as our state variable for monitoring, our inferences are limited to assessing population trends in terms of changes in spatial distribution. However, occupancy tends to be positively correlated with abundance (Gaston et al. 2000), so we expect our results to also provide insights into animal abundance, though this extrapolation requires caution. Occupancy most closely corresponds to abundance when sampling is at the scale of animal home ranges (Linden et al. 2017). In our study, camera spacing in relation to home ranges of muntjac and sambar indicates that animals detected at different sites are probably different individuals. Thus, occupancy estimates should be reliable proxies for abundance for these deer species. Gaur and wild pig roam more widely, and our sampling duration was long enough (2-3 months) to allow them to potentially encounter multiple cameras within their home ranges (though probably <3 cameras, based on known seasonal home range sizes in Thailand; Prayurasiddhi 1997). As a result, some animals detected at nearby sites might be the same individuals. Thus, for gaur and wild pig, the relation between occupancy and abundance is likely to be curvilinear (though still positive) and therefore less sensitive to changes in abundance (Steenweg et al. 2018).

RESULTS
Wild pig was the most commonly detected ungulate species, observed at 51-73% of camera v www.esajournals.org sites each year. Muntjac were detected at 57-66% of sites, sambar at 16-22% of sites, and gaur at 13-19% of sites.
Gaur habitat use was most strongly influenced by proximity to saltlicks, with probability of use declining further from saltlicks (Fig. 3). Two of the five covariate models for gaur failed to converge for 2016, so covariate effects were not assessed for that year. For 2012 and 2014, saltlick had highest AIC weights (55 to 72%) and was the only covariate for which ΔAIC was <2 (Appendix S1: Table S4). Beta coefficients (β) for saltlick were significant each year (Table 1). Other variables had little support.
Sambar occupancy was strongly and consistently influenced by elevation and proximity to saltlicks, with decreasing occupancy further from saltlicks and at higher elevations (Fig. 3). In each year, saltlick and elevation were the sole variables within ΔAIC < 2; AIC weights were >78% each year (Appendix S1: Table S5). Additive models combining these two variables were highly supported each year (model weights >78%), with significant β coefficients. The effect of elevation (β elevation = −1.43 to − 1.535) was 1.3-1.4 times stronger than distance to saltlick (β saltlick = −1.099 to − 1.061). Occupancy probability declined sharply with rising elevation, from 0.44 at 250 m elevation, to 0.02 at 1000 m (Fig. 3).
For wild pig, distance to village, distance to saltlick, and NDVI emerged as the main drivers of habitat use, though not each variable had support in each year. Model uncertainty was particularly high in 2012, with models involving these three variables having comparably low AIC weights (8-18%; Appendix S1:  Fig. 3), but this village effect was not supported in the final year. Saltlick had substantial support in 2012 (sum of AIC model weights containing saltlick = 43%) and 2016 (AIC weight 52%; Appendix S1:  Table S6), with occupancy decreasing further from saltlicks, as with gaur and sambar (Fig. 3). However, the effect of saltlicks on pigs was weaker than for gaur or sambar (Fig. 3). NDVI influenced pig habitat use in 2012, with a positive β coefficient suggesting preference for denser forest over more open habitats.
Muntjac differed from the other ungulates in that the null model, ψ(.), had substantial support each year (ΔAIC 0.77-1.88; Appendix S1: Table S7). None of the variables in the top model set for muntjac were statistically significant (β coefficients crossed 0; Table 1). These results suggest that muntjac distribution was not strongly influenced by any of the predictors we measured.
We obtained 6-18 observations of gaur herds, and 41-65 observations of wild pig herds, each year. Mean herd size of gaur increased slightly over time, from 2.8 individuals per herd (SD 0.83) in 2012, 3.7 (SD 1.7) in 2014, and 4.0 (SD 2.3) in 2016, but this trend was not significant statistically (H = 1.386, df = 2, P = 0.50). Mean herd sizes of wild pig in each year were 17.1 (SD 16.2), 19.1 (21.1), and 13.0 (12.5) and also did not differ statistically (H = 3.014, df = 2, P = 0.222). We examined possible drivers of the increase in wild pig occupancy by comparing models with the colonization parameter as a function of year and ecological or anthropological covariates (ψ(trend), γ(year + covariate), p(effort)). We also compared a no-covariate model: ψ(trend), γ(year), p(effort). The saltlick model did not converge, nor did further additive models with multiple covariates. The most highly supported model was the no-covariate model, with model weight of 36% (Table 2). Although most covariates had some support (ΔAIC < 2), standard errors of β estimates were large, indicating high variability in effects of these variables. Thus, although wild pigs clearly colonized new sites over time, there was no apparent selection for particular types of new sites with regard to the environmental or anthropogenic variables we examined.

DISCUSSION
The ungulate species we studied are among the most frequently poached large mammals in South-East Asia (Harrison et al. 2016). They also constitute main prey of endangered tigers, leopards, and dholes Sunquist 1995, Simcharoen et al. 2018). A critical task accompanying efforts to recover predator-prey communities in South-East Asia is monitoring the status of ungulate populations. Our study identified ecological features that underlie the distribution of ungulate species recovering from overhunting, but also revealed worryingly slow population expansion in most species.

Occupancy patterns
Smaller-bodied ungulate species tend to persist more widely than larger species in the heavily hunted forests of South-East Asia (Vongkhamheng et al. 2013). Our results conformed to this pattern, with muntjac and wild pig occurring at 60-80% of sites but sambar and gaur occupying <25%. Wild pig is particularly resilient to hunting pressure due to its high rate of reproduction (Choquenot et al. 1996). Habitats in the park vary in terms of vegetation type and elevation, but all are suitable for gaur and sambar (Lekagul and McNeely 1988) and it is reasonable to expect the distributions of these species to be much more extensive than we observed. Observations of elder park rangers indicate that sambar and gaur were indeed widespread in the parks before intensive commercial poaching began about 30 yr ago (personal communication with park rangers).  v www.esajournals.org Proximity to saltlicks was a consistently strong influence on the occurrence of wild pig, sambar, and gaur in MWKL (Fig. 3). The effect of distance to saltlick was particularly profound for gaur and sambar: Occupancy rates were just 13.2% and 8.5%, respectively, at sites that were 5 km from a saltlick, and less than 3% at 10 km away (Fig. 3). The effect on wild pigs was less pronounced, with expected occupancy 78% at 5 km and 55% at 10 km (Fig 3). The plant-based diets of herbivores are typically deficient in sodium and other minerals; these minerals are thus obtained from other sources such as saltlicks (Emmons and Stark 1979). The minerals obtained at saltlicks also function to eliminate plant toxins (Kreulen 1985), and provide females with calcium required for gestation and lactation (Kovasc 2005). Wild pigs have broad diets that include occasional animal matter (Oliver and Leus 2008), probably rendering pigs less dependent on saltlicks than sambar or gaur.
Other studies have found saltlicks to influence the distribution and movements of tropical ungulates, including gaur (Prayurasiddhi 1997) and elephants (Chanard et al. 1998). But in contrast to our study, sambar distribution in nearby Huai Kha Khaeng was not influenced by proximity to saltlicks (Simcharoen et al. 2014). There are 37 known saltlicks in MWKL, with a density of about 3 per 100 km 2 . The density of saltlicks (including mineral-rich springs) in Huai Kha Khaeng is about four times higher (12 per 100 km 2 ; data in Prayurasiddhi 1997). The relative scarcity of saltlicks in MWKL might heighten their influence on sambar distribution compared to Huai Kha Khaeng.
Sambar occupancy in most years was associated strongly with low elevations below 1000 m (Fig. 3, Table 1). In 2016, for example, 85% of detections (11/13 camera sites with sambar presence) were below 500 m, and 0% were above 1000 m. The preference of sambar for low elevations was also observed in nearby Huai Kha Khaeng (Simcharoen et al. 2014). Lowland habitats in western Thailand contain important food sources used by sambar that are scarce in montane forests >1000 m, particularly grass and bamboo (Steinmetz et al. 2008).

Population trends and their causes
Of the four ungulate species groups we monitored, only wild pig exhibited a noticeable and statistically significant increase. There are 45 villages within about 5 km from the edge of the parks, and poaching is concentrated closer to villages (unpublished data, MWKL National Parks). Poaching pressure in MWKL is monitored monthly by patrol rangers and is measured in terms of encounter rates with poaching camps, poaching hides, shotgun shells, and snares. Poaching pressure declined by 41% during the study, from a mean of 7.3 encounters/1000 km patrolled in 2012, to 4.3/ 1000 km in 2016 (unpublished data, MWKL National Parks). This might explain why distance to village exerted a negative effect on wild pig occupancy in the initial years of the study, but not in the final year (Table 1). Other studies in the region have also observed significant increases in wild pig occupancy and abundance over similar time periods after poaching declined (Steinmetz et al. 2014). Wild pigs have high reproductive rates relative to the other species in our study, allowing their populations to rebound more rapidly from overhunting (Pepin et al. 2017).
The distributions of muntjac, sambar, and gaur did not increase over the 6-yr study period. These three species are less resilient than wild pig, having slower reproductive rates and lower intrinsic rates of increase. Elsewhere in Thailand, and in Malaysia, sambar have also been particularly slow to recover compared to other ungulate species (Kawanishi et al. 2014, Steinmetz et al. 2014) possibly from excessive poaching of primeage sambar males (targeted for their large antlers) that disrupts the mating system of this species . Our study thus further confirms that sambar do not recover easily from overhunting in South-East Asia. But muntjac and gaur have increased (in both occupancy and abundance) from overhunting within a similar 5-to 6-yr time frame in nearby Thai reserves with similar habitats (Steinmetz et al. 2010(Steinmetz et al. , 2014, so why not at MWKL? One possibility is that our monitoring efforts were insufficient and we failed to detect animals that were there. Although our occasion-specific detection probabilities (p) were indeed somewhat low for all species (<0.28), our extended sampling period (K = 8-12 occasions) resulted in a high probability of detecting each species if present. The probability of detecting a species at v www.esajournals.org least once if present (p*), after K sampling occasions (p* = 1 − (1 − p) K ) ranged from 0.85 to 0.98 each year, with an average across years of 0.92 for muntjac, 0.93 for pig, 0.93 for sambar, and 0.90 for gaur. Thus, our monitoring effort was likely sufficient to detect changes in occupancy. Whereas maximizing detection probability, p, is critical in classic capture-recapture sampling, the focus for occupancy-based monitoring should be achieving high p*, as this increases precision of occupancy models and power to detect trends (Steenweg et al. 2016).
A second possibility is that poaching pressure remained high, perhaps due to lax enforcement efforts by the parks. But this is unlikely, as ranger patrolling effort inside the parks nearly tripled during the course of this study, from 3711 km in 2012 to over 9500 km in 2015 and 2016. Current patrolling effort is actually near the maximum possible, with 11 teams operating >15 d per month. Further, poaching of ungulates does not appear to be severe: between 2012 and 2016, six muntjacs, five pigs, four sambars, and seven gaurs were known to be killed (unpublished data, MWKL National Parks). This amounts to 0.8-1.4 animals of each species killed per year, across an area of almost 1200 km 2 . Although this offtake is likely an underestimate (many poaching incidents might go undiscovered), even if the actual rate were five times higher, it would amount to just 4-7 individuals per species per year. This is likely to be <10% of the extant population of each species, which seems like a low offtake relative to what these ungulate species can likely sustain given their high intrinsic rates of population growth (Hone et al. 2010, Steinmetz et al. 2010. For example, in Europe, red deer (Cervus elaphus), which are similar in body size to sambar, have sustained annual hunting rates of 25% of their populations since the 1980s, yet deer abundance has increased by over 50% (Burbait _ e and Csányi 2010). Further, most poaching in MWKL is focused on small mammals (squirrels, porcupines), junglefowl (Gallus gallus), and reptiles (unpublished data, MWKL National Parks), not ungulates. In sum, the poaching rates of ungulates during our study seem too low to be suppressing their population growth.
A third possibility is that past poaching pressure has induced Allee effects, which inhibit per capita recruitment among remnant ungulates, impeding population growth. Although current levels of ungulate poaching in MWKL are apparently low, this was not the case historically. In the 1970s and 1980s, much of the area was under various resource concessions, including logging, corn cultivation, and Dipterocarp oil extraction. Villages were established to supply and sustain the workers. Roads were built inside the parks to link these concessions to outside towns. These developments facilitated intensive commercial poaching that supplied burgeoning urban markets. Open wildlife markets existed in at least three nearby towns, some operating until 1997 (MWKL rangers, personal communication). The main product sold was game meat-particularly sambar, muntjacs, serow (Capricornis sumatraensis), gaur, wild pig, and primates. Quantities traded were large enough to require dedicated freezers at markets. This trade caused severe declines in poached wildlife in these parks, demonstrated by the complete elimination of banteng (Bos javanicus), and near extirpation of langurs (Trachypithecus phayrei) and gibbons (Hylobates lar).
Small populations of animals can often grow rapidly due to density-dependent effects that enhance survival and reproduction at low density (Fryxell et al. 2014). But if population size drops too low, the per capita growth rate can decline instead (Allee et al. 1949), impeding population recovery. Such Allee effects can be induced, for example, by reduced mating interactions at low density (Courchamp et al. 1999). Gaur and sambar distributions in MWKL could be so sparse that adult females and males rarely meet, thereby reducing reproduction rates and suppressing population growth.
To explore this idea further, we mapped our camera trap records specifically for gaur herds, ignoring single males or bachelor groups. Adult female gaurs occur in herds, so herds are the key reproductive units of a gaur population and underlie its capacity to increase. Gaur herds occurred at just 4-9% of sample locations each year (Fig. 2). Moreover, the few herds that do exist are isolated from each other within MWKL by distance, and from potential connections to gaur populations in adjacent protected areas by high ridges (Fig. 2). Single male gaur range widely over both low and high elevations, but v www.esajournals.org herds tend to stay in the lowlands below 1000 m elevation (Steinmetz et al. 2008).
Making things worse, the size of the gaur herds in MWKL is unusually small, which could further constrain their reproduction. Our camera trap videos of passing gaur herds revealed atypically small herd sizes of 2.8-4.0 (range of averages each year), with a maximum of 7. In contrast, in nearby Huai Kha Khaeng, where a healthy gaur population of over 300 animals occurs, typical herd size was 6-10 animals (Prayurasiddhi 1997). A main function of herding behavior of large ungulates is protection from predators (Fryxell et al. 2007). Dholes, leopards, and tigers prey on gaur, especially the calves (Karanth and Sunquist 1995). If small herd sizes allow exceptionally easy access to gaur calves, then predation from these carnivores could suppress gaur reproduction. Such a predator-mediated Allee effect operating through reduced herd size was implicated in the lack of population recovery of sable antelope (Hippotragus niger) in South Africa (Owen-Smith et al. 2012) and caribou (Rangifer tarandus caribou) in Canada (McLellan et al. 2010). It is possible that such Allee effects are inhibiting population growth of gaur, sambar, and perhaps muntjac, resulting in the flat population trends observed (Fig. 4). Although direct threats from poachers are currently not severe, past poaching has perhaps rendered populations too sparsely distributed to increase.

Management implications
Commercial poaching in South-East Asia is rampant (Gray et al. 2018), and many populations of mammals are in decline (Di Marco et al. 2014). In this context, preventing further declines of wildlife can be considered a partial conservation success. Thus, although occupancy of most species in MWKL did not increase, the stabilization of their distributions is a positive outcome that represents an initial step toward recovery. Nonetheless, we expected population expansion after five years, not just stasis. One lesson from our study is that management efforts to reduce present-day pressure on depressed populations of ungulates might not be entirely fruitful due to historical contingencies, in this case historical poaching. This legacy of poaching might be overriding the effects of current protected area interventions such as patrolling. Learning the environmental history of a project site is thus important for understanding present-day conditions and biological responses to management.
An implication of our study is that intensified law enforcement cannot guarantee recovery of ungulates, especially species such as sambar and gaur that are less resilient to poaching than wild pig. Indeed, gaur, muntjac, and other species have recovered from overhunting elsewhere mainly through interventions other than law enforcement, including collaborative management partnerships with park residents (Steinmetz et al. 2006), and outreach targeting social norms in surrounding villages (Steinmetz et al. 2014). But in those cases, numbers of remnant animals were apparently sufficient to allow recovery once overhunting was alleviated. In contrast, in MWKL, severe overhunting in the past may have pushed populations below a critical threshold, so low levels of current poaching make little difference. Wildlife recovery in MWKL might be constrained further by ecological factors, such as the low density of mineral licks.
We believe more direct interventions such as habitat improvement and population augmentation are needed to increase ungulate distribution and numbers in MWKL, particularly of sambar and gaur. Habitat improvement would entail creating grasslands and artificial saltlicks. The availability of grassy clearings creates favorable conditions for ungulates and supports higher densities than closed forest (Bhattarai and Kindlmann 2012). Saltlicks provide essential nutrients that can improve health and reproductive rates of ungulates. Another benefit of adding saltlicks and grasslands is that, because they act as discrete attractants for surrounding animals, they might facilitate social interactions and increase the probability of mating interactions among remnant animals, thereby helping alleviate one of the Allee effects we suspect is suppressing population growth.
In addition to habitat improvement, we recommend direct augmentation of captive-bred ungulates into the wild. The small population of 9-10 tigers in MWKL has a low reproduction rate, probably due to prey scarcity (WWF-Thailand, unpublished data). Habitat improvement will assist recovery of extant ungulates eventually, v www.esajournals.org but the results will likely be too slow for the immediate requirements of this struggling tiger population. Prey augmentation would more rapidly boost prey distribution and numbers in the short-term. Sambar should be prioritized, as they are a favored prey of tiger (Simcharoen et al. 2018) and are amenable to captive breeding. In the empty forests of South-East Asia, traditional activities such as patrolling and law enforcement may not be sufficient to recover prey populations quickly enough to ensure the viability of tiger populations. More direct interventions like prey augmentation might be required in many parks in Asia that are struggling to conserve their populations of large herbivores and carnivores.