Collaboration for conservation: Assessing countrywide carnivore occupancy dynamics from sparse data

Assessing the distribution and persistence of species across their range is a crucial component of wildlife conservation. It demands data at adequate spatial scales and over extended periods of time, which may only be obtained through collaborative efforts, and the development of methods that integrate heterogeneous datasets. We aimed to combine existing data on large carnivores to evaluate population dynamics and improve knowledge on their distribution nationwide.

| 919 VAN DER WEYDE Et Al. Utilizing data from many different sources, however, requires analytical tools able to integrate heterogeneous datasets, a common problem in ecology.
In recent years, several analytical methods have been developed to assess species occupancy, persistence and colonization.
Occupancy, the probability that a species is present at a particular site (MacKenzie et al., 2002), is widely used where abundance or density is more difficult to estimate, and it provides the ability to measure species-, site-or survey-specific detection probabilities (Broms et al., 2016). A growing body of occupancy studies has demonstrated the potential to develop multi-species approaches (Broms et al., 2016), to consolidate data from multiple survey methods to improve parameter estimates (Linden et al., 2018;Nichols et al., 2008) and to use single-visitation data to estimate detection parameters (Lele et al., 2012;Peach et al., 2017). Yet, where data are sparse in space and time, and heterogeneous in structure, there remains a need to incorporate these aspects into a single analytical framework. This creates opportunities to harness existing datasets from a variety of sources that could be used to improve species management by assessing the impact of various environmental and anthropogenic factors on species distributions and dynamics at large spatial scales and over extended periods of time.
Anthropogenic and environmental factors such as human distribution and activities, land use, habitat characteristics, prey availability, rainfall, as well as intra-guild competition and predation have been shown to affect carnivore distribution, abundance, movement and survival (Qi et al., 2020;Ripple et al., 2014). For instance, while the larger, conspicuous lion (Panthera leo) and spotted hyaena (Crocuta crocuta) are often confined to protected areas (Bauer et al., 2015;Mills & Hofer, 2013), the smaller and more mobile cheetah (Acinonyx jubatus), African wild dog (Lycaon pictus), leopard (P. pardus) and brown hyaena (Parahyaena brunnea) demonstrate an ability to also occupy agro-pastoral regions (Boast & Houser, 2012;Lindsey et al., 2004). These regions vary in both human and livestock densities, thus influencing the presence of large carnivores through variations in persecution risk, prey populations and water availability. The influence of anthropogenic and environmental factors on carnivores may act locally or at larger spatial scales and effects may manifest seasonally or over longer periods of time.
Particularly in arid and seasonally fluctuating environments, spatiotemporal changes in prey availability and surface water can have important consequences on movement and distribution of carnivores.
Although many wildlife species are resilient to small fluctuations in anthropogenic and environmental factors (Qi et al., 2020), large carnivore range and resource requirements render them particularly vulnerable to extinction where drastic habitat changes occur (Ripple et al., 2014).
Few sub-Saharan countries still support an intact large carnivore guild with viable populations. Across most of southern Africa, habitats and populations have already been fragmented severely (e.g. Riggio et al., 2012). Botswana is an exception, with strong carnivore populations that occur contiguously across most of the country (Ray et al., 2005;Winterbach et al., 2014). Recent studies have highlighted Botswana's value in ensuring connectivity of lion populations across neighbouring countries, an essential requirement for long-term genetic viability (Cushman et al., 2018). Similarly, Botswana provides a critical linkage between cheetah populations scattered across the southern African sub-region (Weise et al., 2017). It follows that maintaining within-country connectivity will be paramount for both national and international conservation efforts. There is a need, therefore, to improve carnivore distribution knowledge in Botswana, particularly beyond the well-known populations residing in protected areas (PAs), and to identify key areas that are important for their long-term survival, while maintaining connectivity between populations.
Our objective was to formulate a dynamic occupancy model that could integrate data from a range of survey methods and species to investigate spatial and temporal factors influencing large carnivore distribution in Botswana. We predicted that major drivers like human population, livestock density, rainfall and land use would be key determinants affecting occupancy of all species, while seasonal changes in primary productivity would mainly influence short-term persistence and colonization. Additionally, we aimed to use our model to predict the probability of occurrence of each species countrywide, thus also covering areas that lack empirical data.

| Study area
Botswana covers approximately 582,000 km 2 and has a growing population of more than 2.3 million people (Statistics Botswana, 2014). Located on a plateau, its mean altitude above sea level is approximately 1,000 m (Burgess et al., 2006). Botswana is relatively flat and, except in the north with its seasonally inundated Okavango Delta-Chobe River ecosystem, it is a mostly arid country dominated by savanna biomes and the Kalahari Desert. Several large salt pans (Makgadikgadi-Ntwetwe-Sua ecosystem) that span over 20,000 km 2 occur in the centre of the country. Rainfall is highly seasonal with a dry season between May and October and a wet season between November and April (Batisani & Yarnal, 2010). The vegetation is determined by a strong north (>600 mm) to south (<200 mm) rainfall gradient (Cooke, 1985).
Nearly 40% of Botswana comprises formally protected areas like National Parks and Game Reserves, as well as Wildlife Management Areas (WMAs). While the latter are not gazetted as formally protected by the government (Twyman, 2001), wildlife-oriented activities (e.g. photographic and hunting tourism) are the main form of land use in WMAs (Parry & Campbell, 1990). Much of the remainder of the country is used for agro-pastoral activities (primarily cattle or crop farming) on communal, leasehold or freehold land (National Resource Services, 2002). Botswana's urban population is mostly concentrated in large cities, towns or villages, in the eastern part of the country, while the rural population is spread among remote villages, farms and cattle posts (Statistics Botswana, 2014).

| Data collection and processing
We gathered presence-only and presence/absence data for six large carnivores (African wild dog, brown hyaena, cheetah, leopard, lion and spotted hyaena) from 31 datasets collected by 19 independent researchers or organizations over the period May 2010-April 2016.
We considered four methods of data collection: camera trapping, spoor surveys along transects, telemetry data from very high frequency (VHF) and global positioning system collars (GPS) collars, and opportunistic sightings. Where telemetry data originated from translocated individuals, we only used locations after 1 year from translocation to ensure that the data represented settled individuals.
We divided Botswana into a 10 × 10 km grid in ArcGIS (ESRI,v10.3), resulting in a total of 5,985 grid cells. We asked each researcher/organization to provide data on any or all of the six large carnivores detected during their research in any of the grid cells: for each grid cell, species were reported as either detected or not detected with a specific method ( Figure 1). We pooled the data into 6month periods coinciding with the wet and dry season, respectively, providing a total of 12 sampling sessions for the entire study. We included in the analysis grid cells that were surveyed with any of the four methods during at least one of the 12 sessions.
We calculated sampling effort in each grid cell during each session for camera trap and spoor surveys to account for the differences in effort between datasets and cells. For camera traps, we summed the total number of recording days or camera days across all active devices per grid cell. Sampling effort varied between 0 and 108 days/session. For spoor surveys, we summed the total sampling distance across all transects surveyed in a grid cell during a particular sampling session. Sampling effort varied between 0 and 100.45 km/ session. As telemetry records and opportunistic sightings were available as presence-only data, we could not estimate sampling effort for these methods.

| Predictors of large carnivore occurrence
We collated data on five covariates: land use, human population size, livestock density, rainfall and NDVI that we considered to have strong justification for influencing carnivore distribution (Henschel et al., 2016;Midlane et al., 2014;Mkonyi et al., 2018;Rich, Davis, et al., 2017;Winterbach et al., 2014). We used ArcMap (ESRI,v10.3) for data compilation, digital mapping and spatial analysis. Layers were in Africa Albers Equal Area Conic projection with D WGS 1984 datum. Maps of each layer are provided in the Figure S1.
We pooled land use types across Botswana into four categories, which best represented the majority of the country. Land use types were as follows: formally protected national parks and game F I G U R E 1 Location of grid cells (10 × 10 km) that were surveyed across Botswana and those showing detections of at least one of the six large carnivores. Land use divided into four categories: protected areas (PA), wildlife management areas (WMA), commercial farms (CMF) and communal grazing areas (CGA) are also shown distributed across the country reserves used for recreational purposes (PAs); wildlife management areas used for recreational and some agricultural activities (WMAs); freehold and leasehold commercial farms and ranches commonly used for agro-pastoral activities (CMFs); and communal grazing, arable land and settlements (CGAs). We used the majority tool in zonal statistics in ArcMap to identify the dominant land use for each cell.
To calculate human population size, we used data from the 2011 National Census, as this was the most recent information available.
We created a human population index map by summing population size values of all the settlements falling into each grid cell. To account for the fact that human influence may extend beyond a particular grid cell, values for single cells were the sum of its value plus the value of its adjoining eight cells. Although this will potentially overestimate population size, it provided a more realistic layer for showing relative variation in human population density across the country. We checked our map against the Global Population Count and Human Footprint maps (CIESIN, 2016) to ensure our covariate layer was suitable.
We calculated livestock density from aerial surveys conducted by the Department of Wildlife and National Parks (DWNP) of Botswana in 2012, 2013 and 2015. We included cattle, donkey, sheep, goat and horse data. As aerial surveys did not cover the same spatial extent each year, we pooled data from the three surveys to create a single livestock density map covering the entire country. Livestock survey data, recorded as the number of livestock per 1 × 1 km raster pixels, were first averaged for overlapping grids and mean values were then created for each of our 10 × 10 km grid cells for the analysis.
Rainfall data for each season were available from rainfall stations (47-60 per session) located across the country. We interpolated these data using inverse distance weighting (IDW), resulting in a 1 × 1 km raster file. We subsequently rescaled raster pixels to our 10 × 10 km grid cells. We then combined seasonal rainfall data into a single layer that represented the average rainfall over the entire study period. This represents long-term climatic conditions that we assumed affects the large-scale distribution of all species and allowed us to assess the broader impact of above-ground water availability on carnivore occupancy along the north to south rainfall gradient that characterizes Botswana.
We utilized the normalized difference vegetation index (NDVI), which represents seasonal-specific changes in vegetation productivity, as a proxy for herbivore occurrence (Mueller et al., 2007;Pettorelli et al., 2009). We deemed this "vegetation greenness index," which results from a combination of factors including rainfall, vegetation type, soil type and ambient temperature as appropriate to represent seasonal changes in herbivore distribution and availability. For each year, we obtained NDVI data using MODIS13A3 imagery from EarthData obtained from http://lpdaac.usgs.gov. We used 1 × 1 km resolution MODIS/NDVI data (Didan, 2015) from June and December as representative months for the dry and wet season, respectively. We averaged and rescaled the 1 × 1 km pixel data to match our 10 × 10 km grid cells. We standardized all continuous covariates to a mean of zero and a standard deviation of one, and we checked for correlation between covariates using variance inflation factors.
We limited our modelling of covariates for specific parameters for biological reasons and limitations due to our sparse data (Table 1). Human factors (e.g. population size and livestock density), land use and rainfall are likely to affect large-scale distributions of species across the country (Mkonyi et al., 2018;Weise et al., 2017).
Similarly, rainfall can drive net primary productivity and has previously shown to affect occupancy (Henschel et al., 2016); however, mean rainfall should not directly affect extinction or colonization probabilities. Rather these parameters are largely driven by human factors or seasonal variation in productivity, as measured by NDVI (Fernandez et al., 2016).

| Modelling framework
For our dataset, we only had summarized data of whether a species was detected or not, together with the total survey effort. Lele et al. (2012) showed that occupancy and detection probability can be accurately estimated from single survey data, provided there is at least one continuous detection covariate that is not also used to estimate occupancy. Peach et al. (2017) extended this approach to a dynamic occupancy model.
Given that our data spanned 6 years, with two seasons in each year for a total of 12 sessions, and that we were interested in modelling occupancy trends over that period, we used a dynamic occupancy model to explicitly model local persistence and colonization (MacKenzie et al., 2003;Royle & Kéry, 2007). We extended the model by Peach et al. (2017) in two ways. First, we took advantage of the fact that we had data from four survey methods (although not necessarily all of them at the same site) and therefore used a multimethods model (Nichols et al., 2008), which is analogue to a multiple independent-observer model (Alldredge et al., 2006). Second, we formulated the model as a multi-species model to jointly model the six species together (Dorazio & Royle, 2005) and implemented it in a Bayesian framework.
Our response variable consisted of presence/absence data y ijtm for each of the six species i = 1, 2, … I and survey grid cell j = 1, 2, … J, collected during each of the twelve sampling periods t = 1, 2, … T with any of the four sampling method m = 1, 2, … M.
We modelled the true occupancy state z ijt (z = 1: occupied, z = 0: unoccupied) over time t, where ψ ijt = Pr(z ijt = 1) is the occupancy probability. We further defined the probability that a site remains occupied from time t to t + 1 ("persistence") as ɸ ijt and the probability of colonization for an unoccupied site as γ ijt . We modelled the true occupancy state z ijt as a Markovian process based on the initial state z ij1 : We modelled parameters as a function of covariates using a dot product of two vectors: Xocc j , Xper jt and Xcol jt are covariate vectors for site j at time t and βocc i , βper i and βcol i are vectors of speciesspecific regression coefficients with occ = occupancy, per = persistence and col = colonization: For our multi-species model, we modelled each regression coefficient hierarchically as random factors e.g.
where μ is the mean and σ 2 is the variance of a particular coefficient β across all species (Kéry & Royle, 2008).
For the two survey methods for which survey effort was quantifiable (i.e. camera traps and spoor surveys), we modelled the detection probability p ijtm as a power function based on the detection probability per unit effort r ijm and the total sampling effort k jtm for method m at grid cell j during sampling period t (Peach et al., 2017).
For survey methods that returned presence-only data and for which survey effort was not quantifiable (i.e. telemetry and opportunistic sightings), we set r ijm = 1 and k jtm = 1 (resulting in p ijtm = 1) for grid cells where the species was detected and k jtm = 0 (resulting in p ijtm = 0) for grid cells where the species was not detected.
The observed data y ijtm is then the outcome of a Bernoulli process: We implemented the model in JAGS (Plummer, 2003) through software R (R Development Core Team, 2019). We computed three chains with 80,000 iterations, a burn-in of 40,000 iterations and a thinning rate of 40. We assessed convergence visually and considered the Gelman-Rubin statistic with a value of <1.1 acceptable (Gelman & Rubin, 1992). Code for running the model in JAGs is provided as Figure S2.
We mapped changes in occupancy over time by subtracting the predicted occupancy values for sampled grids of the last session from predicted occupancy values of the first session. Lastly, we predicted occupancy nationwide, that is also for unsampled grid cells using the same dynamic model for unconditional estimates of psi. Lions and African wild dogs were the most and least detected species, respectively ( Table 2).

| RE SULTS
Occupancy was negatively correlated with livestock density for all species except the brown hyaena (Table 3). This species had strong evidence for higher occupancy in commercial farms than PAs, whereas African wild dog, leopard and lion had lower occupancy.
For all species, occupancy probability did not differ between protected areas, WMAs or communal grazing areas. A strong negative correlation with human population size was evident only for lion, indicating minimal opportunity for lions to persist in densely humanpopulated areas. Rainfall positively correlated with occupancy of the two largest species, lion and spotted hyaena, while it was negatively correlated with brown hyaena occupancy (Table 3), highlighting the adaptation of the latter to arid environments.
The relationship between persistence and colonization and the three covariates varied considerably across the focal species (Table 3). While persistence of brown hyaena and spotted hyaena was positively associated with human population size, the probability of colonizing new areas by these species had a negative relationship with human population. In contrast, leopard persistence was TA B L E 1 Summary of covariates considered to influence large carnivore occupancy (ψ), persistence (φ), colonization (λ) and detection probability (p)  We used the results from our model to generate nationwide predictions for each species.

| D ISCUSS I ON
Our study is one of the first attempts to model carnivore occupancy on a national scale. The resulting knowledge is crucial for improving management efforts, particularly for a country like Botswana that holds some of the largest remaining carnivore populations in Africa (Weise et al., 2017;Winterbach et al., 2014) and thus plays a critical role in their conservation. We extended the single-visit multi-season occupancy model of Peach et al. (2017) by adding a multi-species and multi-method component. This allowed us to utilize existing data from different sources and generated with different survey methods.
This novel approach was necessary to increase and expand sampling resolution and cover several ecosystems, thus improving our ability to identify factors influencing species distribution across years and over large spatial scales. A similar analysis would not have been possible using data collected by single researchers who generally focus on smaller study areas spanning a few thousand square kilometres (Cozzi et al., 2013;Van der Weyde et al., 2018;Weise et al., 2018).
We show that occupancy was mainly explained by livestock density and land use regime. Except for brown hyaena, occupancy was lowest in commercial rangelands and where livestock density was highest. One explanation is that, except for brown hyaenas that scavenge almost exclusively (Maude, 2005), our focal species regularly prey upon livestock and, therefore, suffer from retaliatory persecution (Boast et al., 2015;Weise et al., 2018) and may generally prefer to avoid human-dominated areas. Brown hyaena may thus take advantage of reduced intra-guild competition (Kent & Hill, 2013;Winterbach et al., 2017). A similar scenario has been hypothesized for cheetahs (Durant, 2000;Wachter et al., 2011; but see Scantlebury et al., 2014;Swanson et al., 2014), which are known to persist on farmlands (Boast & Houser, 2012;Weise et al., 2017).
Our results indicated a strongly negative relationship between human population size and lion occupancy. TA B L E 2 Number of grid cells with detections during at least one survey for each of the six focal carnivore species from four survey methods. Detections in cells between species were not mutually exclusive found to occupy PAs and WMAs. This highlights the need to safeguard these conservation areas and to enhance their connectivity to maintain genetically healthy lion populations (Cushman et al., 2018).
Additionally, we found that communal rangelands are also utilized by many large carnivores and may serve as important habitats. Leopard occupancy was predicted to be high across the country, and this likely reflects its highly adaptable nature and tolerance of humanaltered environments (Stein & Hayssen, 2013). Our grid size could also have influenced our findings, as truly unsuitable patches (e.g. city centres) may indeed be smaller than those used in this analysis.
Due to the ecological plasticity that leopard exhibit, other predictors such as degree of anthropogenic persecution or intra-guild competition may be more useful for predicting distribution. Finally, we observed a positive relationship between rainfall and occupancy of lions and spotted hyaenas, the two largest species. This may be best explained by their preference for abundant, large prey such as Cape buffalo (Syncerus caffer caffer), zebra (Equus quagga) and wildebeest (Connochaetes taurinus), which commonly occur in large numbers in wetter, highly productive environments (Selebatso et al., 2018).
We found little consistency in the factors explaining carnivore persistence in already occupied areas and colonization of new areas.
The negative effect of increasing human density on leopard persistence is not surprising considering that leopards are commonly reported to prey on livestock and likely to suffer from persecution in Botswana (Schiess-Meier et al., 2007;Selebatso et al., 2008), nor is the positive persistence of brown hyaena, a species known to cope well in human-dominated landscapes (Kent & Hill, 2013).
However, the positive correlation between human density and persistence of spotted hyaenas was unexpected, and we have no clear explanations for this. Spotted hyaenas are highly adaptable and able to adjust their behaviour to benefit from humans (Kolowski & Holekamp, 2008;Yirga et al., 2016). In Botswana, however, they are generally perceived as a conflict species and not tolerated close to human habitation or livestock dwellings. As expected, high human density correlated negatively with all large carnivores' ability to recolonize an area. The effect of NDVI on persistence was more variable among species. Spotted hyaena had a higher probability of persistence in more productive environments, in agreement with our findings that this species occurred mostly in the northern regions. By contrast, lion and cheetah appeared to persist better in drier areas.
This may be due to the highly seasonal rainfall that occurs in the northern areas of the country where prey abundance and distribution are seasonally variable (Bennitt et al., 2014;Kotze et al., 2020) forcing these species to adapt their movements in response to seasonal changes. We acknowledge that other factors such as population density of each species, which can lead to different degrees of intra-and inter-specific interactions (Durant, 2000;, might be as (or even more) important than environmental and anthropogenic factors to understand patterns of persistence and colonization. Although we did not attempt to account for cooccurrence patterns mainly due to limited sample sizes, we encourage doing so in future occupancy applications.
Our model predicted substantial portions of the country are suitable to maintain viable populations of the six large carnivores.
Occupancy results indicated a distinct spatial separation between the two hyaena species, even though both species exhibited colonization into each other's mainstay range. Climate change-induced aridity may enable the dryland-dwelling brown hyaena to expand its distribution northward, while the creation of artificial watering holes in Botswana's arid south might facilitate spotted hyaena range expansion, as it already has with other wildlife (Selebatso et al., 2018).
We, however, lack empirical information to formally test this hypoth- Combining a large number of datasets allowed us to produce the first estimates of carnivore distribution dynamics across Botswana, but the sparse nature of the data also resulted in some limitations.
Dynamic occupancy models require repeated sampling of the same site to estimate probabilities of persistence and colonization (MacKenzie et al., 2003). There are few long-term carnivore research projects in Botswana, and, where they exist, they are often limited to small spatial scales, or areas that are more accessible or preferred by researchers (e.g. wildlife versus human-dominated areas), potentially biasing our findings. Consequently, the 6-year dataset lacked temporal replicates for some areas and species. Furthermore, there were few data for eastern Botswana, as well as the country's south-east, far south and north-west. This limits our inferences and increases uncertainty, as reflected by the model's large credible intervals. The temporal resolution of covariate data also potentially restricted our ability to detect significant changes. For example, we were only able to include a single estimate of human density and livestock abundance across the entire study period, although we expect that these changed during the study, leading to stronger impacts on wildlife distribution. Our use of NDVI as a proxy for prey availability may also be limited in identifying the strong influence that prey is known to exercise on carnivore distribution (Rich, Davis, et al., 2017;Winterbach et al., 2013). To accurately portray temporal trends and identify major factors influencing carnivore occupancy, strategic long-term research in select locations would improve estimate precision. Even so, the need for constructive collaboration remains, through data sharing and otherwise, to be able to produce national population estimates for single or multiple species, while also identifying major threats or positive conservation measures. Our work highlights the necessity, and the feasibility, of collaborating on large-scale assessments of wildlife populations.
Only rarely will any single organization or government agency manage to collect an adequate amount of data allowing assessment of wildlife distribution and population size at large scales.
Even if data exist, the propensity of not sharing information, heterogeneous data structures, low consistency in repeat surveys and unequal sampling effort can present considerable challenges (Hampton et al., 2013;Rich, Davis, et al., 2017;Weise et al., 2017).
Developing species distribution and population models that integrate heterogeneous data while accounting for imperfect detection are advancing rapidly (Linden et al., 2018;Peach et al., 2017), and our approach may assist in addressing the challenges of heterogeneous data common in wildlife ecology. Although we show that accommodating such data heterogeneity is possible, we stress the importance of regular coordination and collaboration among researchers and organizations to achieve the best possible conformity in data collection. Additionally, we highlight the value of maintaining heterogeneous land use types and livestock distribution in supporting carnivore conservation. Here, we focused on large carnivores, but our approach can be extended to any other species or taxa with data available at extensive temporal and spatial scales (e.g. Chase et al., 2018;Keeping, 2014). We stress that nationwide approaches and collaborations such as the one presented in this study are crucial to drawing inferences at adequate biological and spatial scales, to develop research priorities, and for the implementation of evidence-based conservation and management strategies.

ACK N OWLED G EM ENTS
The

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interests arising from this study.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data associated with this paper are available in the Dryad Digital Repository https://doi.org/10.5061/dryad.ghx3f fbp9.