Getting the big picture: Landscape‐scale occupancy patterns of two Annamite endemics among multiple protected areas

The Annamite mountains of Vietnam and Laos are a global biodiversity hotspot harboring several threatened endemic species. Conservation efforts to protect these endemics are hampered by a lack of knowledge on their ecology and distribution. We conducted landscape scale camera‐trapping across six study areas in the Annamites to assess distribution patterns of two endemics: the Annamite dark muntjac species complex Muntiacus rooseveltorum/ truongsonensis and the Annamite striped rabbit Nesolagus timminsi. We used a Bayesian single‐species occupancy framework to estimate occupancy as a function of ecological and anthropogenic factors. Our study showed that Annamite dark muntjac was predominantly found at higher elevations (>1000 m) and in areas that were more inaccessible to people and had lower surrounding village density. Annamite striped rabbit exhibited both positive and negative responses to elevation among study areas, with no clear response to the anthropogenic covariates. Our results showed that covariate responses varied among the study areas when random effects were included on study areas. We discuss the application of random effects to investigate species occupancies across large spatial scales, and the risk of not accounting for variation among study areas. Our prediction maps provide the first comprehensive overview of the distribution of these endemic species across a substantial part of their range and can be used to help stakeholders focus conservation efforts on priority areas.

Bayesian single-species occupancy framework to estimate occupancy as a function of ecological and anthropogenic factors. Our study showed that Annamite dark muntjac was predominantly found at higher elevations (>1000 m) and in areas that were more inaccessible to people and had lower surrounding village density. Annamite striped rabbit exhibited both positive and negative responses to elevation among study areas, with no clear response to the anthropogenic covariates. Our results showed that covariate responses varied among the study areas when random effects were included on study areas. We discuss the application of random effects to investigate species occupancies across large spatial scales, and the risk of not accounting for variation among study areas. Our prediction maps provide the first comprehensive overview of the distribution of these endemic species across a substantial part of their range and can be used to help stakeholders focus conservation efforts on priority areas.

K E Y W O R D S
Annamite dark muntjac, Annamite striped rabbit, camera-trapping, endemic species, Laos, occupancy modeling, tropical rainforest, Vietnam

| INTRODUCTION
The Annamite mountain range of Vietnam and Laos is a global biodiversity hotspot harboring exceptionally high levels of mammalian endemism (Baltzer et al., 2001). Several Annamite endemic mammals were only recently discovered by science (Giao et al., 1998;Do et al., 1994;Jenkins et al., 2005;Van Dung et al., 1993), highlighting how little is known about the biodiversity in this region. Unfortunately, all mammals are threatened by the widespread unsustainable hunting pressure that is prevalent across the region (Brook et al., 2014;Gray et al., 2018;Harrison et al., 2016). Hunting in the Annamites is primarily accomplished by the setting of indiscriminate wire snares (Belecky & Gray, 2020;Gray et al., 2018), and over the past several decades, this threat has decimated populations of terrestrial mammals in the region. For instance, both the saola Pseudoryx nghetinhensis and the large-antlered muntjac Muntiacus vuquangensis are listed as Critically Endangered on The IUCN Red List for Threatened Species, and are rapidly approaching global extinction (Timmins et al., 2016(Timmins et al., , 2020. To protect remaining populations of these endemics, it is imperative to reduce snaring in priority areas. However, the ability of conservation stakeholders to effectively utilize limited resources to reduce snaring pressure in key areas has been hindered by our poor understanding of the ecology and distribution of these species. Filling this knowledge gap remains an important step toward implementing more effective conservation actions in the region.
Among the Annamite endemics, the Annamite striped rabbit Nesolagus timminsi and Annamite dark muntjac species complex Muntiacus rooseveltorum/ truongsonensis are among the most understudied. The Annamite striped rabbit is a forest-dwelling lagomorph that was only discovered by science in the mid-1990s (Averianov et al., 2000;Surridge et al., 1999). It is currently listed as Endangered on The IUCN Red List for Threatened Species as a result of substantial snaring-driven declines that are believed to have occurred across its range (Tilker et al., 2019). The little work that has been conducted on the Annamite striped rabbit indicates that hunting pressure is a primary driver of its distribution in some areas, and that populations are strongly depressed in key strongholds . The Annamite dark muntjac complex is believed to contain two, or perhaps more, muntjac species (Timmins & Duckworth, 2016a, 2016b. Because of its unresolved taxonomy, the dark muntjac complex is currently listed as Data Deficient, though there is little question that its populations have been depressed across the Annamites through unsustainable hunting. The importance of accounting for imperfect detection has become a fundamental tenet of ecological studies in recent decades Kéry & Schmidt, 2008;Lahoz-Monfort et al., 2014). Through the use of hierarchical occupancy models, researchers are able to explicitly model the detection process, thus providing a more accurate estimate of species occurrence (MacKenzie et al., 2002;Kéry & Andrew Royle, 2016). Not surprisingly, this analytical advancement is particularly important for surveys involving rare or cryptic species, and has become standard practice for cameratrapping studies involving elusive rainforest mammals (Johnson et al., 2009;Sollmann et al., 2017). Assessing the distribution of tropical mammals across large spatial scales often requires joint analyses of data from multiple independent study areas, thus necessitating strategies for data pooling (Gelman & Hill, 2006). One option is to estimate independent parameters for each study area (nopooling). The no-pooling approach, however, is often not possible due to sparse data, particularly for rare species. Furthermore, differences in covariate ranges among the studies might lead to a biased or limited understanding of species-covariate relationships along covariate gradients and within study areas. In contrast, a complete-pooling strategy combines all data without differentiating among study areas, thus assuming that the ecological relationships between the species and the covariates remain constant across all study areas. While this approach may help provide parameter estimates, especially for rare species with few detections, it might be inappropriate in some contexts. For example, the complete-pooling strategy does not account for ecological differences among study areas and variation in species responses to environmental factors, which is often the main interest of ecological studies (Kéry & Schmidt, 2008;Tobler et al., 2015;Rich et al., 2017;Johnson et al., 2020). Ignoring the local relationships between species occurrence and covariates may therefore bias occupancy estimates and predicted habitat associations. An increasingly popular compromise is partial-pooling, in which the effects at individual study areas are considered to stem from a shared underlying distribution. Due to this shared underlying distribution, partial pooling allows inferences for species with sparse data for which heterogeneity between study areas (or other levels) is expected. Partial pooling is implemented in models via random effects (Gelman & Hill, 2006).
Here, we used landscape camera-trapping data and an occupancy modeling framework to estimate study area-level occurrence for Annamite dark muntjac and Annamite striped rabbit across six study areas in the Annamites. To explore how species responses to covariates differ among study areas, we used two modeling approaches: First, we used fixed effect models in which species responses to covariates remained constant across all study areas. Second, we allowed speciescovariate relationships to vary among study areas by including a study area random effect on the covariate responses and intercepts. The first approach reflects the common strategy of complete-pooling of data from different areas to increase detections, but at the cost of ignoring potential variation among study areas. The second approach utilizes a partial-pooling strategy, allowing heterogeneity within the species-covariate relationships among study areas while permitting areas with insufficient data to "borrow" information from those with more data (Korner-Nievergelt et al., 2015). We then used results from the study area-specific random effect approach to predict species occurrence in our six study areas, thus providing the first landscape-scale overview on the distribution of these little known Annamite endemics.

| Study area
We carried out systematic landscape-scale cameratrapping surveys in six study areas located within the northern and central Annamites. In the northern Annamites, we surveyed Pu Mat National Park (NP), Vietnam. In the central Annamites, we surveyed four areas in Vietnam: Bach Ma NP, Phong Dien Nature Reserve (NR), the contiguous Hue and Quang Nam Saola NRs, and Song Thanh NP. We also surveyed one area in Laos including the eastern section of Xe Sap National Protected Area (NPA) and an ungazetted forest to the south of Xe Sap NPA near the village of Ban Palé ( Figure 1). In our analyses, we considered the two Saola NRs and the Xe Sap NPA/Palé sites as single study areas, respectively (Table 1).
All areas are characterized by steep and rugged terrain with mean slopes ranging from 19 to 27 and elevations ranging from 200 to 1800 m asl (Baltzer et al., 2001;Le et al., 2001;SFNC, 2000). Annual precipitation ranges from 1800 to 8000 mm and the predominant habitat is evergreen broadleaf forest with a multitiered closed canopy (Baltzer et al., 2001;Long, 2005;SFNC, 2000). All areas are managed by local and national government authorities, with the exception of the Ban Palé area, which does not have protected area status.
There are considerable differences among the areas in their accessibility and levels of human use. The Vietnamese study areas are surrounded by local communities living within the reserve buffer zones. Pu Mat NP and Song Thanh NP have lower surrounding village densities than the other three Vietnamese study areas, though there are five villages located within the southern part of Pu Mat NP core zone. In Phong Dien NR, there are four hydropower dams within the protected area, and a major road running through the park that provides relatively easy access to central and southern parts of the protected area. The Saola NRs and Bach Ma NP are also bisected by major roads. In contrast, the Laos study areas are more remote and characterized by lower surrounding population densities.
Among the six study areas, Pu Mat NP and the Saola NRs are the only two areas with moderate to high levels of snare-removal. Since 2011, WWF-Vietnam and provincial authorities have maintained professional ranger teams within the Saola NRs, with the primary goal of finding and removing wire snares . Over the past 10 years, these teams have removed more than 110,000 individual snares from the reserves . Pu Mat NP has also received support for snare removal activities from a collaboration between the park and local nongovernmental organizations (NGOs), though intensive patrolling has been implemented over a shorter time span. Bach Ma NP, Phong Dien NR, and Song Thanh NP have not benefited from the high levels of external technical and financial investment that characterize the other study areas, and as a result, on-the-ground enforcement efforts in these areas have been less intensive. In Laos, the eastern Xe NPA/Palé area has received some patrol effort, but the intensity and coverage have not been as high as in the neighboring Saola NRs.

| Systematic camera-trapping survey
We conducted systematic camera-trapping surveys between November 2014 and March 2019 (Table 1). All surveys followed a standardized approach, with stations spaced approximately 2.5 km apart across the survey areas and operational for approximately 60 days (see Abrams et al., 2018 for more details on the study design). At each station, two nonfacing cameras (Reconyx 850 or Panthera V6 white-flash units) were set within a 20 Â 20 m plot. Detections from the two cameras were treated as from a single station in our modeling framework. Cameras were placed on trees 20-40 cm above the ground and were programmed to take a three-photo burst with no delay. We set up a total of 368 camera-trap stations across all six study areas (Table 1). The resulting camera trap data were managed using the R package camtrapR (Niedballa et al., 2016).
To prepare the camera-trap data for occupancy analyses, we produced detection histories for Annamite dark muntjac and Annamite striped rabbit using an occasion length of 5 days. We assumed spatial independence between stations based on available information on species home range sizes. Although there is little information on the home range size of the Annamite dark muntjac, studies on other muntjac species have indicated a home range of approximately 1 km 2 (McCullough et al., 2000;Odden & Wegge, 2007). For Annamite striped rabbit, a recent study in the Saola NRs showed that the movement patterns of individual rabbit are likely below 500 m . Thus, we assumed the average home range sizes for both species to be smaller than the minimum distance between our cameratrap stations. We nevertheless acknowledge that the occupancy probabilities we present may be more appropriately interpreted as probabilities of habitat use due to the point sampling characteristics of camera traps and uncertainties about species' home range sizes and densities (Efford & Dawson, 2012). Steenweg et al. (2018), F I G U R E 1 Map of study areas however, showed that occupancy estimates and inferences from point sampling are robust and irrespective of spatial grain of analysis irrespective of species ecology. Because the sampling period (73 ± 15 days) was small in comparison to the life span of our target species, we also assumed that the true occupancy state did not change during the survey.

| Model covariates
We included covariates on both detection and occupancy probabilities. For detection probability, we used a measure of survey effort, defined as the number of active camera-trap nights per station and occasion. This was performed to account for failure of individual camera trap units and incomplete occasions.
We modeled occupancy as a function of three covariates: elevation, a least-cost path measure of accessibility, and village density surrounding the study areas. We considered elevation to be a primarily environmental covariate because previous studies have shown fine-scale differences in habitat structure along elevational gradients in the Annamites (Long, 2005). We considered the least-cost path and village density covariates to be proxies for anthropogenic pressure in our study areas.
We extracted elevation values at camera-trap stations from a Shuttle Radar Topography Mission (SRTM) 30-meter digital elevation model (Farr et al., 2007) using the raster package in R (Hijmans, 2020). Studies have shown Annamite dark muntjac occurrence to be positively associated with elevation (Long, 2005;Timmins & Duckworth, 2016a, 2016b. Less is known about the elevational ranges of the Annamite striped rabbit, though there is some evidence that the species occurs across a range of elevations in different areas (Dang et al., 2001;. We created a least-cost path layer for each study area using Tobler's hiking function (Tobler, 1993) implemented in the R package movecost (Alberti, 2019). With this approach, we calculated the travel time from major roads to each of the sampling stations. The output is a measure of time needed to reach a given point in the study areas, which we interpret as a measure of remoteness. To capture historical levels of accessibility, we conducted the analyses using road layers from three different time periods spanning the past 20 years. The initial temporal starting point was taken from a road layer from the year that the sampling occurred (see the survey period in Table 1); we then used road layers from the previous 10 and 20 years for the subsequent layers. The final leastcost path layer was created by averaging the accessibility for the three time periods (see  for more information on the least-cost path calculations). The least-cost path raster was produced at a 200-m resolution, and sampling station values were extracted using the raster package (Hijmans, 2020).
We created village density layers to serve as a proxy for human pressure around the protected area. Previous studies, both in the Annamites  and other tropical regions (Harrison et al., 2016), indicate that hunting pressure is positively associated with village proximity. For this covariate, we first created a point layer representing the villages surrounding each protected area. Village locations were derived from a combination of government sources and Google Earth satellite images. All point layers were then combined into one layer from which we generated a village density heatmap raster using the heatmap plugin in QGIS v3.14 (QGIS Development Team, 2020). Because our field experience indicated that all camera-trap stations to be subject to some level of hunting pressure, we chose a radius of 20 km for the village heatmap so that the resulting raster layer covered all sampling stations. The village density layer was also produced at the 200-m resolution. We plotted the values of covariates at camera-trap stations to visualize the overlap of covariates ranges among study areas ( Figure S4, Supporting Information). Generally, the elevation and least-cost path value ranges were similar among study areas, while the village density range slightly varied. All covariates were scaled to have a mean of 0 and standard deviation of 1 to allow for direct comparison of model estimates and improved parameter estimation. A Spearman rank correlation showed that no two covariates were highly correlated (|r| < 0.44; Figure S5). To ensure consistency of spatial scales for the occupancy predictions, we resampled the 30-m elevation raster to the same 200-m resolution so that it was at the same spatial scale as the least-cost path and village density layers.

| Modeling framework
To model species occurrence, we used single-species occupancy models (MacKenzie et al., 2006) implemented in a Bayesian framework. We fitted three occupancy models for each species: (i) a fixed effect model which included three covariates as fixed effects without accounting for study areas, (ii) a random effect model which included a random effect of study area on each of three covariate and intercept (without fixed effects), and (iii) a null model which did not include covariates on occupancy. Model structure was identical between species to facilitate comparisons between the two species (see Data S1, Methods).
Models were run using the R package ubms version 1.0.2 (Kellner, 2021). We used the default vague priors for all models (normal distribution with mean = 0 and SD = 10). We ran four parallel Markov chains with 5000 iterations each, and discarded the first 2500 iterations as burn-in. Because ubms uses the Stan modeling languageand therefore, Hamiltonian Monte Carlo sampling instead of Gibbs sampling-relatively low numbers of iterations are sufficient to achieve model convergence and stable parameter estimates (Carpenter et al., 2017). Model convergence was assessed by checking that the Rhat statistics equal 1 and by visually examining the model traceplots (Gelman et al., 2014). To determine the significance of the covariate effects, we generated 95% and 75% credible intervals of the posterior distribution. We considered a coefficient to have strong support if its 95% Bayesian credible interval (BCI) did not overlap zero, and moderate support if the 75% BCI did not overlap zero.
We compared and ranked the three different candidate occupancy models of each species in ubms using expected log pointwise predictive density (elpd) as a measure of predictive accuracy of the models (more specifically, elpd assesses how well a model predicts future data, assuming that the future data come from the same distribution as the observed data). Elpd was calculated using leave-one-out cross-validation for pairwise model comparisons (Vehtari et al., 2017). To assess model support relative to the top model, we also calculated the pairwise differences in elpd (Δelpd) between models and the top model and the standard errors thereof (SE [Δelpd]). We considered a model to have less support than the top model if its absolute difference in elpd was greater than the standard error of that difference (Doll & Jacquemin, 2019;Kellner, 2021).
To compare results from the fixed and random effect models, we generated covariate response curves (marginal effect plots) for Annamite dark muntjac and Annamite striped rabbit, and then predicted occupancy for each study area under each modeling approach. We therefore produced two occupancy prediction maps for each study area for both Annamite dark muntjac and Annamite striped rabbit. We then calculated the difference in predicted occupancy probabilities between both models at all raster cells to highlight the discrepancies between the spatial predictions under each approach.
Finally, we generated estimates of percentage of area occupied (PAO) for each model and study area by calculating the occupancy probabilities at each raster cell for each posterior sample, converting these to binary occupancy status at each cell as the outcome of random binomial trials, and aggregating these to PAO for each study area and each posterior sample. By repeating this process for all posterior samples, we obtained distributions of estimated PAO for all study areas.

| Camera-trapping effort and species detections
In total, we obtained data from 366/368 camera-trap stations totaling 43,452 camera-trap nights (Table 1). Two stations did not provide data because the cameras failed at one station and the units were stolen at another. Overall, we obtained 110 independent detections of Annamite dark muntjac from 41 stations and 173 independent detections of Annamite striped rabbit from 41 stations (Table 1).

| Annamite dark muntjac
Model selection showed that the random effect model had the highest predictive accuracy (measured by elpd), followed by the fixed effect model and null model, respectively. Pairwise jΔelpdj of both the fixed effect model and the null model relative to the random effect model was higher than their SE(Δelpd) (Table S3). Thus, the random effect model was considered to be the most supported model.
In the fixed effect model, we found that elevation had a strong positive effect, the least-cost path had a moderate positive effect, and village density had a strong negative effect on Annamite dark muntjac occupancy (Figure 2a,b). In the random effect model, the association between the covariates and occupancy varied among study areas. Elevation positively affected occupancy in four of the six study areas (Pu Mat NP, Xe Sap NPA/Palé, Saola NRs, and Song Thanh NP), though the significance and strength of the relationship varied (Figure 2a,b). The least-cost path covariate had a strong positive effect in one study area (Pu Mat NP) and a moderate effect in another (Xe Sap NPA/Palé). Overall, village density had less impact on occupancy than the two other covariates, with a moderate negative impact in one study area (Phong Dien NR).
Among study areas, the mean estimated PAO of the Annamite dark muntjac was similar between the fixed and random effect models. The mean estimated PAO ranged from 0.07-0.29 based on the results from the fixed effect model, and 0.06-0.31 based on the random effect model. Both models showed that the Xe Sap NPA/Palé area had the highest mean predicted PAO and Bach Ma NP had the lowest (Table S4). Differences in the predicted occupancy on a pixel level between the fixed and random effect models ranged from À0.48 to 0.28 ( Figure 3). Overall, the two models had similar results in estimated occupancy probability for Annamite dark muntjac across the six study areas. However, the random effect model showed more spatially pronounced predicted distribution patterns, with areas having both higher and lower predicted species occupancies, whereas in the fixed effect model the variation within study areas was less pronounced (Figures 3 and S2).

| Annamite striped rabbit
In the model selection, we found that the random effect model had the highest predictive accuracy (measured by F I G U R E 2 Standardized beta coefficients and response curves. The beta coefficients (a, c) indicate the mean and 95% Bayesian credible interval (BCI) as a thin line, 75% BCI as a thick line showing covariate effects from the fixed effect and random effect models on occupancy of both study species. Red bars indicate strong effect where the 95% BCI does not overlap zero, black bars indicate moderate effect where the 75% interval does not overlap zero but the 95% interval overlaps zero, and gray bars show effects in which the 75% BCI overlaps zero (no effect of covariate). The response curves (b, d) showing the response of species occupancy probabilities relative to the original covariate values (bold lines are posterior means, ribbons are 95% BCIs). The black dashed line indicates the mean of responses across all study areas elpd), followed by the null and fixed effect models, respectively. Pairwise Δelpd of the null model was smaller than its SE(Δelpd) ( Table S3) and was therefore considered to have similar model support as the random effect model. The fixed effect model showed that elevation had a moderate positive effect and least-cost path had a moderate negative effect on Annamite striped rabbit occupancy. We found no support for the effect of village density on species' occupancy (Figure 2c,d). In the random effect model, elevation was the only covariate that affected the Annamite striped rabbit occupancy. The effect of elevation varied among the study areas, with a strong positive effect in Song Thanh NR, a moderate F I G U R E 3 Predicted occupancy probability for the Annamite dark muntjac at the six study areas (in rows) from the fixed effect model (fixed effect column) and the random effect model (random effect column). Disparity column shows the difference between predicted occupancy for the Annamite dark muntjac from the fixed and random effect models. PAO column shows the percentage of area occupied (PAO) estimation for the Annamite dark muntjac positive effect in the Saola NRs, and a strong negative effect in Pu Mat NP (Figure 2c,d).
The mean predicted PAO of the Annamite striped rabbit based on the fixed effect model was similar among the six study areas (0.13-0.16; Table S4). However, the random effect model estimates of PAO showed considerable differences. The random effect model predicted that the Saola NRs had the highest mean PAO (0.26 ± 0.07) and Pu Mat NP the lowest PAO (0.05 ± 0.03) ( Table S4). As with the Annamite dark muntjac, the occupancy predictions for the Annamite striped rabbit were more spatially pronounced in the random effect model than in the fixed effect model (Figures 4 and S2).
F I G U R E 4 Predicted occupancy probability for the Annamite striped rabbit at the six study areas (in rows) from the fixed effect model (fixed effect column) and the random effect model (random effect column). Disparity column shows the difference between predicted occupancy for the Annamite striped rabbit from the fixed and random effect models. PAO column shows the percentage of area occupied (PAO) estimation for the Annamite striped rabbit 4 | DISCUSSION

| Importance of including study areaspecific covariate effects
For both Annamite dark muntjac and Annamite striped rabbit, the random effect model had higher model support than the fixed effect model-despite being a more complex model with more parameters-indicating that it explained the data better (Table S3). The fixed effect model showed an overall moderate positive association of Annamite striped rabbit occupancy with elevation. This positive association primarily comes from pooling data across all study areas in the fixed effect model, and is misleading in some areas; for example, in Pu Mat NP, Annamite striped rabbit were only recorded at lower elevations. Within this context, predications made in these areas using pooled data may not present an accurate assessment of Annamite striped rabbit occurrence with respect to elevation. The random effect model, on the other hand, was able to disentangle the different responses of Annamite striped rabbit to elevation among the different study areas, thus accounting for heterogeneous habitat associations for the species. For example, we found positive (Song Thanh NP and Saola NRs) and negative (Pu Mat NP) relationships of Annamite striped rabbit occupancy to elevation in different study areas. The predicted distributions using the random effect model were therefore based on more supported habitat associations and showed differences in the mean PAO among the study areas (Table S4).
We also found that the random effect model provided a better fit than the fixed effect model for Annamite dark muntjac. However, the difference between the outputs of the random effect model and the fixed effect model of the muntjacs was less pronounced in comparison to the Annamite striped rabbit. Nevertheless, the random effect model showed variation in the Annamite dark muntjac's study area-specific relationship with the covariates. For example, although the fixed effect model showed a strong positive effect of elevation on the Annamite dark muntjac occupancy pattern in all study areas, we were able to describe study area-specific nuances in this response using the random effect model. For instance, in Pu Mat NP the response of Annamite dark muntjac to elevation was better described as a moderate positive effect of elevation (random effect model) than a strong positive effect (fixed effect model). Therefore, even when species relationships to landscape factors are broadly similar among study areas, allowing for study area-specific heterogeneity in the strength of these relationships may provide a more accurate picture of occurrence.
In this study, we included random effects on all covariates and the intercept, which may be unnecessary in studies if there is no a priori reason to assume random variation in the responses to all covariates across the study areas. In situations where the covariate responses are similar among study areas, accounting for minor differences via random effects or fixed/interaction effects may not outweigh the disadvantages associated with increased model complexity. Should there be heterogeneous responses to certain covariates only, it would be possible to specify models with fixed effects of some covariates and random effects on others. To identify which approach is more suitable, one could compare parameter estimates from single area occupancy models first, and if they are found to differ strongly, account for this variation by including a random effect in a combined model. Alternatively, one can formulate different combined models (corresponding to competing hypotheses about habitat associations at different study areas) and compare them via model selection techniques (Aho et al., 2014;Johnson & Omland, 2004).
Another consideration is that random effect models often show higher uncertainty in parameter estimates compared to fixed effect models (Kéry & Schaub, 2012). These higher uncertainties are a consequence of the added complexity and additional model parameters associated with random effect models compared to fixed effect models. Researchers should be mindful of the implications of random effect models and adjust their data analysis strategies accordingly (Bolker et al., 2009).
Overall, our results suggest that complete pooling of data and failing to account for heterogeneous responses to covariates at different study areas may lead to failure to identify meaningful, study area-specific covariate relationships. This, in turn, may potentially lead to erroneous conclusions about species distributions and thus misguided conservation decisions. Our results indicate that researchers should be cautious when combining data across multiple study areas in a complete-pooling approach since it assumes that species-habitat relationships are consistent across these areas. Ultimately, when and how researchers utilize random effect models to account for heterogeneity in species responses will be context-specific and will depend on the model species, the sampling areas, the data quality and quantity, and the overall aims of the study.
We discuss the implications of our results on species ecology and conservation below. Because the random effect models fit our data better than the fixed effect models, we focus only on the study area-specific random effect model output.

| Annamite dark muntjac
Elevation was the most important driver of Annamite dark muntjac occurrence in our study areas. The fact that elevation had a significant positive effect on Annamite dark muntjac occupancy in four out of the six study areas generally supports earlier findings that suggest that it occurs at higher elevations (Long, 2005;Timmins & Duckworth, 2016a, 2016b. However, our results also show that the overall effect of elevation on Annamite dark muntjac occurrence is complex. Clearly, the Annamite dark muntjac tends to be more common at higher elevations, but the importance of elevation as a driver of its distribution varies across its range. Explaining these elevational patterns remains challenging, in part because general bioclimatic patterns in the Annamites are complex and still poorly understood, and in part because elevation is often conflated with anthropogenic pressures. From a purely ecological perspective, our results provide support for the idea put forward by Duckworth (2016a, 2016b) that the Annamite dark muntjac is predominantly found in forests that are more mesic, which tend to be found in higher elevations. The clearest indication for this finding is shown in Song Thanh NP, where wet evergreen forest is conspicuously confined to high elevation areas along the border with Laos (Long, 2005); all of our Annamite dark muntjac records from Song Thanh NP were from this area. Interestingly, we did not find a strong effect of elevation in either Phong Dien NR or Bach Ma NP. One possible explanation for this finding is that forest toward the more mesic end of the Annamitic evergreen forest spectrum extends to lower elevations in these areas. One possible explanation for Phong Dien NR is that, in contrast to most other areas in the Annamites, the high elevation areas are easily accessible from surrounding local villages, whereas the less accessible lower and mid elevation in the north-west areas where we detected dark muntjac (ranging from 174 to 709 m a.s.l.) are more inaccessible ( Figure S3); studies in other regions have shown that mammals may persist in suboptimal habitat where threats are lower (Kerley et al., 2012).
We were surprised that the least-cost path covariate did not show a stronger response among our study areas. Interestingly, the only areas where least-cost path had a strong effect were the two largest areas that we surveyed (Pu Mat NP and Xe Sap NPA/Palé). The least-cost path values for these areas indicate that these study areas have levels of remoteness that far exceed the other study areas ( Figure S3). Our observations in the field support this finding: It was only in Pu Mat NP and Xe Sap NPA/Palé that our teams noticed a declining trend in human signs and snares in the more remote areas. It is therefore possible that we only detected a strong signal for the least-cost path covariate in the two study areas which contain sufficiently remote regions.
The only study area where village density had an effect was Phong Dien NR, which showed a moderate negative relationship between the covariate and Annamite dark muntjac occurrence. The result from Phong Dien NR is broadly consistent with our expectations that population pressure would be inversely correlated to hunting pressure. We are unable to explain why village density did not have a more pronounced effect on Annamite dark muntjac distribution in our other study areas. It is possible that village density by itself is too coarse to serve as a meaningful proxy for hunting pressure in some areas within the Annamites. For example, the village density covariate does not capture potential cultural differences that may influence hunting behavior among different ethnic minority groups (Harrison et al., 2016). It is also possible that actions on the ground have a greater influence on current hunting pressure than village density by itself: the Saola NRs have high village density values, but these areas have what is probably the most intensive snare-removal operation in all of Vietnam .

| Annamite striped rabbit
Elevation was the only covariate with at least moderate support for an effect to Annamite striped rabbit occupancy. It is interesting that the strength and direction of the response to elevation varied among three of our study areas (Figure 2c,d). As with the Annamite dark muntjac, these elevational patterns may be related to the complex bioclimatic conditions in the Annamites. The fact that the Annamite striped rabbit was found in highland areas of Song Thanh NR is consistent with the idea that the species is a wet evergreen forest specialist (Tilker et al., 2019). Perhaps the most unexpected finding was the negative effect of elevation on Annamite striped rabbit occurrence in Pu Mat NP. Previous studies have documented the species at low elevations in Pu Mat NP (Abramov et al., 2016;SFNC, 2000), though to our knowledge, this is the first study to show that the species is either absent or occurs at low occupancies at higher-elevation forest areas within the park. Higher elevation areas of Pu Mat NP are likely to be cooler during the winter months than areas in the central Annamites (SFNC, 2000;Tran et al., 2008), and it may be that these conditions are suboptimal for the Annamite striped rabbit. Future studies are needed to investigate the potential effects of seasonality on Annamite striped rabbit elevational ranges in the Annamites.
We were surprised to find that neither the least-cost path nor village density were informative for Annamite striped rabbit in any of our six study areas. There could be multiple reasons that we failed to find a response with these covariates. It might be that the Annamite striped rabbit is not as sensitive to hunting pressure as larger species such as ungulates, or as noted above, these proxies may not adequately capture the potentially complex spatial patterns of hunting pressure within the Annamites. Different snaring techniques might have different impacts on Annamite striped rabbit populations, for example, if brush snare-fences are substantially more detrimental than single hoop snares; such complexities would not be captured in our least-cost path or village density covariates. It may be that a more sophisticated proxy will need to be developed to fully capture these anthropogenic pressures for the Annamites, including potentially subtle study area-specific differences.

| Conservation of two Annamite endemics
Our predicted species distributions provide the first largescale overview of the potential distribution of these two little-known Annamite endemics. Our results indicate that the distributions of both species are likely spatially restricted within each of the six study areas. The heterogeneity in the responses suggests different drivers of species distributions are at play in the different study areas. Such a spatial restriction likely reflects complex biogeographical patterns as well as the fact that both endemics have experienced population declines. From an applied perspective, the distribution maps can help local management authorities to target their enforcement activities to certain areas within the study areas. Although it would be ideal for entire protected areas to receive sufficient levels of snare removal, in reality, most of the protected areas in our study are far from having the resources, technical capacity, and political support to implement such an approach (with the notable exception of the Saola NRs, where over 110,000 snares were removed the last 10 years; . Therefore, it may be necessary to target snare removal activities to areas of the highest biodiversity value and conservation concern. Moreover, the predicted distribution maps may help to identify other areas within our study areas where these rare species occur. Such information can guide future surveys that seek to gather additional information on these species, for example through the collection of eDNA to assess population genetics (Nguyen et al., 2021).
Interestingly, the highest occupancies for Annamite striped rabbit were in the Saola NRs, where snare removal efforts have been the most intensive and sustained among our six study areas. In fact, the Hue and Quang Nam provincial border areas had the highest occupancy probability for both species. An earlier community occupancy analysis also supported the importance of this area for threatened mammals . Future surveys in the Saola NRs are needed to provide information on the impact of intensive snare removal efforts on recovery of highly depleted populations. We recommend that upcoming patrolling efforts in other protected areas should follow standardized protocols, which were set by WWF Vietnam for their work in the Saola NRs, with all results documented using the Spatial Monitoring and Reporting Tool (SMART). With standardized protocols, it may be possible to identify hunting hotspots, assess hunting trends over time, and evaluate the effectiveness of snare removal activities . Applying this established system in other protected areas would also allow comparisons of hunting pressure among protected areas, potentially allowing for the identification of broader hunting trends, and thus help to set conservation priorities and measure the success of conservation interventions.

ACKNOWLEDGMENTS
Funding for the surveys was provided by the CarBi project, Green Annamite project, German Federal Ministry of Education and Research (BMBF FKZ: 01LN1301A), Leibniz-IZW, Point Defiance Zoo and Aquarium, Safari Club International, Critical Ecosystem Partnership Fund, the National Geographic Society's Committee for Research and Exploration, Ocean Park Conservation Foundation Hong Kong, Mohamed bin Zayed Species Conservation Fund, and GreaterGood. Thanh V. Nguyen received financial support through a Russell E. Train Education for Nature Program fellowship, the Leibniz-IZW, and the Vietnamese Ministry of Science and Technology Program 562. The authors declare no competing interests.

CONFLICT OF INTEREST
The authors declare that there is no conflict of interest that could be perceived as prejudicing the impartiality of the research reported.

AUTHOR CONTRIBUTIONS
Andreas Wilting and Andrew Tilker conceived the ideas. Thanh V. Nguyen, An Nguyen, Andrew Tilker, and Anh C. Dao collected the data. Thanh V. Nguyen and Jürgen Niedballa analyzed the data. Benjamin M. Rawson, Anh Q. H. Nguyen, Trung T. Cao, and Oliver R. Wearn assisted the field surveys. Thanh V. Nguyen, Andrew Tilker, Andreas Wilting, and Jürgen Niedballa led the writing. All authors commented on and reviewed the manuscript.

DATA AVAILABILITY STATEMENT
Underlying data, R codes that support this study are available on the GitHub repository (https://github.com/ EcoDynIZW/Nguyen_2021_CSP.git). Any requests for a full data for academic purpose, please contact t. nguyen@izw-berlin.de or wilting@izw-berlin.de.