Study area

Our study area was the northern section of the Northern Territory of Australia, commonly known as the ‘Top End’ (Fig 1). The region is subject to a wet-dry tropical climate where annual rainfall ranges between ~1400 mm in the north to 500 mm in the south, of which about 90% falls during a pronounced wet season (typically November-May). Tropical savanna extends as far south as the 500 mm annual rainfall isohyet [20], and contains large areas of eucalypt woodlands, eucalypt open forest, and smaller areas of monsoon rainforest, Allosyncarpia forest, sandstone heath, floodplain, wetland and riparian communities [21].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Site locations sampled between 2011 and 2016 using multiple methods.

Location of 333 monitoring sites (black dots) in eight protected areas (1-Gregory National Park, 2-Nitmiluk National Park, 3-Kakadu National Park, 4-Litchfield National Park, 5-Garig Gunak Barlu National Park, 6-Warddeken IPA, 7-Djelk IPA, 8-Fish River Station) in the Top End of the Northern Territory, Australia, from which presence-absence data for mammals, reptiles and birds was collected.

https://doi.org/10.1371/journal.pone.0203304.g001

Data collation

We utilised presence-absence data from the most recent round of a long-running vertebrate monitoring program within Kakadu, Litchfield and Nitmiluk national parks [11], combined with data from surveys in other areas across the study area. In total, we collated repeat detection/non-detection data from 333 sites across five national parks (Kakadu, Litchfield, Nitmiluk, Garig Gunak Barlu, and Gregory), two Indigenous Protected Areas (IPA; Warddeken and Djelk), and one privately owned conservation reserve (Fish River Station) (Fig 1). Sites within each reserve stratify all major vegetation types within dryland habitats of the Top End. We limited our dataset to surveys conducted between 2011 and 2016 to provide an up-to-date snapshot of species occupancy.

All sites were visited once and surveyed for mammals, birds and reptiles using a standard protocol [8], although the timing of visits to sites within a year varied between parks. Sites consisted of a 50 x 50 m quadrat (hereafter referred to as the trapping quadrat) containing: 16 Elliott traps placed equidistantly around the perimeter; eight cage traps placed at the corners and mid-way along each side of the quadrat; and three 20-L pitfall traps with 10 m of drift fence. Cage and Elliott traps were baited with a mixture of peanut butter, honey and oats, set at dusk and checked and closed at dawn. Pitfall traps were checked for reptiles and mammals at dawn, midday and dusk. Cage, Elliott and pit trapping was conducted for three or four days/nights at each site. Surveys were completed under Charles Darwin University Animal Ethics Permit A13026, and all sampling procedures were specifically approved as part of obtaining the field permit.

At a subset of sites, diurnal bird surveys (225 sites), nocturnal searches (228 sites), and camera trapping (163 sites) were conducted across a 100 x 100 m quadrat centred on the trapping quadrat (S1 Table). Three ten-minute bird surveys were conducted per day over three days. Additionally, one fifteen-minute nocturnal search by spotlight was conducted on each of three nights after dark targeting reptiles, mammals and birds, without using call-playback methods to increase detection rates. Five Reconyx® motion-sensor cameras were deployed during the four-day survey following the methods of Gillespie et al. [22] and remained operational for five weeks to detect mammals. Camera images were processed by a minimum of two experienced personnel. We acknowledge these sampling methods are inadequate for amphibians and bats. Due to the expectation of infrequent records for these poorly sampled animal groups we excluded them from the analysis.

Detection histories

We generated detection histories for species observed at least once during the study: a 1 represented the detection of a species, a 0 represented a non-detection. We collapsed detection histories so that the sampling occasion was one day/night for live trapping, bird surveys and nocturnal searches. This resulted in three to four repeat sampling occasions to sites for these methods. For camera trapping, we collapsed detection histories so that the sampling occasion was one week, resulting in five repeat visits to sites with this method. Separate detection histories were generated for each sampling method. For example, a detection history of [0,1,0,0,1] with camera trapping and [1,0,0,0] with cage trapping meant a species was detected (one or more individuals) in the second and fifth week with cameras, and only in the first day/night with cages.

Environmental predictor variables

We compiled a list of covariates that were: 1) thought to broadly influence occupancy and detectability of birds, mammals and reptiles across the study region; and, 2) available at the broad spatial scale of the study. These included climatic variables (mean annual rainfall, maximum temperature in the warmest quarter, minimum temperature in the coldest quarter), terrain variables (elevation, terrain ruggedness, distance to perennial creeks), environmental variables (soil clay content, vegetation cover), and fire variables. Fire variables included point-based metrics (number of fires, time since fire) and neighbourhood-based metrics (fire extent, proportion burnt and fire patchiness) to capture the spatial homogeneity of fires across space and time [12]. We acknowledge that other covariates such as predator density or habitat structure may also drive occupancy of mammals, reptiles and birds. However, such factors could not be included in our modelling because where data exists it is patchy, hence could not be mapped for the study region. We obtained raster layers of the covariates and re-sampled each one in ArcGIS software at 100 m resolution (S1 Fig) before extracting their value at sampling sites. For a detailed description of the covariates and justification for their selection see Table 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Summary of the 13 climatic, topographic and fire covariates considered in occupancy and detectability modelling.

https://doi.org/10.1371/journal.pone.0203304.t001

Data exploration

We tested site covariates for collinearity (S2 Table) and discarded one of a pair if the Spearman correlation coefficient was greater than 0.7 [29]. As a result, we discarded minimum temperature in the coldest period (correlated with rainfall) and proportion burnt (correlated with fire frequency) from any further analysis.

Occupancy and detectability modelling

We modelled occupancy and detectability for each species assuming the result from any given survey was the outcome of two binomial processes acting simultaneously [30, 31]: 1) the probability a species was present at a site (ψ) over long time periods; and, 2) the probability a species was present within the site and observed in any given survey visit (ρ). The first binomial process, occupancy, was considered a Bernoulli random variable, such that: (1) where ψ is the probability of species occupancy at site i and z_i denotes the true state of occurrence. We then assumed the observation of a species at site i is the outcome of another Bernoulli random variable, with the product of z_i and detection probability ρ_ij the success probability. The detection/non-detection data y_ij observed at site i during survey j can be described as: (2) where y_ij is the observed ‘presence-absence’ data and ρ_ij is the probability of detection. We modelled the influence of covariates on occupancy and detectability using the logit link function in Eqs 1 and 2.

Model fitting

Prior to model fitting, we identified seven species (Arnhem Land Rock Rat, Black Wallaroo, Gehyra pamela, Pseudothecadactylus lindneri, Oedura gemmata, Banded Fruit-dove, Chestnut-quilled Rock-pigeon) that have ranges restricted to the Arnhem Plateau, which is only a small fraction of the study region. We created a ‘plateau’ layer as a proxy for this geographic range by placing a 5 km buffer around the Arnhem escarpment, and only used data collected from within this area when fitting models. This meant that our estimates of occupancy were relative to the expected geographic range of species rather than to all sites in the study region, with predictions of the seven species above confined to our plateau layer.

We fitted occupancy and detectability models using the unmarked package [31] in R [32]. We assumed that occupancy could depend on any combination of environmental predictor variables, while detectability was modelled as a function of the sampling method (for those species with multiple detection methods), time since fire, fire frequency, and terrain ruggedness (Table 1). Our rationale for including these variables in the detectability component of the model is that we expected fire to influence vegetation structure and composition at sites and therefore possibly effect the ability of observers to sight individuals. Terrain ruggedness might also influence how well individuals can be sighted during surveys by limiting the line of sight of the observer.

We appended the detection histories for each species and modelled the sampling method in unmarked as an observational-level factor. We assumed a closed population over the survey period or that movement of species with large home ranges to and from sites was random [33]. We allowed for quadratic terms in all covariates in the occupancy model except for fire extent and fire patchiness, but did not expect any quadratic terms in the detectability component. All covariates were centred and standardised prior to model fitting [34].

We fitted the most complex model to each species and considered all possible combinations of covariates using the pdredge function in the R package MuMIn [35]. To avoid over-fitting, we constrained the number of covariates in each model at a ratio of 1:10 with the number of sites where a species was detected [36]. We ranked all models by their AIC value [37] and retained only those with a ΔAIC of less than six [38].

Model selection

For each species, we screened our competing set of top models for ecological realism. Reviewing and scrutinizing models so that predictions conform to expert opinion or beliefs is generally encouraged in the species distribution modelling literature [39]. To aid this validation process we generated occupancy maps by predicting the top models to covariate raster layers mapped at 1 km resolution (S1 Fig). We generated 1 km maps from our original 100 m raster layers by reclassifying each covariate layer. Re-sizing rasters was done simply to reduce the image size of each map (<20 MB) for ease of use.

Occupancy maps were screened and discarded if: 1) there was perceived to be a severe mismatch between predicted occupancy and current beliefs about where the species occurs; or, 2) if predicted occupancy unrealistically approached 0 and/or 1 (i.e., boundary estimates), which can result from too few detections [40]. Consideration was given to species with very low occupancy and/or detectability (<0.1) as such estimates should be treated with caution [30]. Models with the highest AIC were accepted as the best model unless occupancy maps failed the validation process, in which case the next-best ranked model was considered.

We assessed model fit of the top-ranked model with sum of squared errors, Freeman-Tukey Chi-squared test, and Pearson’s Chi-squared test [41]. Goodness-of-fit tests were conducted with 1000 simulations. If the top-ranked model for a species failed any goodness-of-fit test it was discarded for the next-best ranked model. If we could not fit a model containing one or more covariates, we modelled constant occupancy across the study area with the sampling method as an observational-level factor (for those species with multiple detection methods). For each species, we calculated the mean occupancy and detectability for each method by averaging across estimates at sampling sites.

Mapping species occupancy and richness

For species containing covariates in their ‘best-model’, we generated occupancy maps across the study region by predicting occupancy to mapped covariate raster layers at 1 km resolution (S1 Fig). Rather than predict occupancy to the entire study area, we limited our predictions to the mainland because offshore islands are disjunct locations where species might respond differently to the environmental covariates. We were also cautious about predicting occupancy well-outside the range of the sampled environmental domain. After comparing the environmental domain at monitoring sites (i.e., range of predictor covariate values) with that of the study area we masked out regions with much higher maximum temperature and fire patchiness, and much lower annual rainfall and vegetation cover, than what was sampled. This restricted our occupancy predictions in some parts of the study region, particularly in the south.

We built relative species richness maps by stacking occupancy maps for those birds, mammals, and reptiles for which we could fit models, and summed the predicted occupancy value in each cell across species [42]. We did not include species with constant occupancy (i.e., null models, or no environmental predictor variables in the model) and note that our species richness maps do not represent true species richness for the study area, but rather the richness of the species we could fit models to.

Occupancy

A total of 242 native species (147 birds, 69 reptiles and 26 mammals) were recorded during surveys of eight protected areas (S3–S5 Tables), representing 66% of the total number of bird, mammal and reptile species sampled over 20 years of the Three Parks Program. We fitted occupancy-detection models to 65% of these, of which a further 9% were discarded because they were considered implausible (S6–S9 Tables). This resulted in occupancy-detection models for 136 species, including 83 birds, 33 reptiles, and 20 mammals. Birds exhibited the most variability in mean occupancy across sampling sites, ranging from 0.02 ± 0.01 (Arafura Fantail) to 0.81 ± 0.08 (White-bellied Cuckoo-shrike) (Fig 2; S6 Table), while reptiles ranged from 0.01 ± 0.01 (Ctenotus storri) to 0.49 ± 0.07 (Carlia amax) (Fig 2; S7 Table). By comparison, the range in mean occupancy for mammals was lower than birds or reptiles, ranging from 0.02 ± 0.01 (Northern Quoll) to 0.30 ± 0.06 (Agile Wallaby) (Fig 2; S8 Table).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Occupancy versus detectability for mammals, reptiles and birds in Northern Australia.

Occupancy (x-axis) versus detectability (y-axis) for mammals, reptiles and birds ± SE (grey bars) estimated with repeat detection/non-detection data collected at 333 sites in Northern Australia between 2011–2016. Different shapes represent the different sampling methods. Detectability for live trapping and spotlighting is estimated over the period of one day/night, while estimates of camera trapping is over the period of one week.

https://doi.org/10.1371/journal.pone.0203304.g002

For birds, elevation was retained in 34% of occupancy models, followed by maximum temperature (27%), fire frequency (22%), fire patchiness (16%), distance to watercourse (16%), fire extent (13%), and time since fire (13%, Fig 3). For reptiles, maximum temperature was retained in 35% of occupancy models, followed by ruggedness (32%), clay content (25%), and fire frequency (21%, Fig 3). For mammals, fire extent was the covariate most commonly included in occupancy models (33%), followed by fire frequency and vegetation cover (both 25%), ruggedness and elevation (both 18%, Fig 3). Of the 136 modelled species, 121 included one or more covariates (S2 Fig). Maps of predicted occupancy and the response of occupancy to covariates are shown for a selection of species in Fig 4, and occupancy maps for all species are presented in S3–S5 Figs.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Environmental covariates commonly retained in occupancy and detectability models.

The proportion of times (y-axis) each covariate (x-axis) was included in occupancy and detectability models for birds (top row), reptiles (middle row) and mammals (bottom row). Site covariates are grouped by landscape, climate or fire. Dark shading indicates the proportion of species that responded negatively to covariates, white shading indicates the proportion that responded positively.

https://doi.org/10.1371/journal.pone.0203304.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Relationship of covariates with occupancy for a selection of species.

Example of three maps of predicted occupancy for Common Rock Rat (top left) Spangled Drongo (top centre) and Carlia amax (top right) across the Top End of Australia. The relationship between the probability of occupancy and selected covariates (scaled) for these three species are presented in the bottom row. Dotted lines represent 95% confidence intervals.

https://doi.org/10.1371/journal.pone.0203304.g004

Detectability

Mean detectability during a day/night of live trapping (mammals and reptiles) or day/night of searches (birds, mammals and reptiles) was lowest for mammals (0.13 ± 0.05), and higher for birds (0.23 ± 0.07) and reptiles (0.28 ± 0.10), although estimates varied widely between species and between sampling methods (Fig 2; S6–S8 Tables). For example, daily detectability of birds (three surveys) ranged from 0.01 ± 0.03 (Collared Sparrowhawk) to 0.74 ± 0.01 (Rufous-banded Honeyeater, S6 Table). Nightly detectability of birds using spotlighting (one survey) was lower, ranging from 0.03 ± 0.02 (Barking Owl) to 0.44 ± 0.21 (Spotted Nightjar). For reptiles, detectability per day/night using the array of pitfall traps ranged from 0.02 ± 0.02 (Heteronotia planiceps) to 0.69 ± 0.11 (Ctenotus piankai; S7 Table), whereas detectability during a night of spotlighting (one survey) ranged from 0.10 ± 0.05 (Gehyra nana) to 0.57 ± 0.22 (Oedura marmorata). Of three reptiles detected using both methods, two were more detectable by spotlighting (S7 Table).

Detection probability for mammals using a night of pitfall trapping was very low, ranging from 0.01 ± 0.01 (Grassland Melomys) to 0.09 ± 0.06 (Red-cheeked Dunnart) (Fig 2; S8 Table). In contrast, nightly detectability using spotlighting varied widely from 0.03 ± 0.03 (Sandstone Antechinus) to 0.67 ± 0.16 (Sugar Glider). Similarly, nightly detectability using the array of Elliott traps ranged widely, from 0.01 ± 0.01 (Delicate Mouse) to 0.57 ± 0.08 (Arnhem Land Rock Rat). The Arnhem Land Rock Rat had the lowest detection probability using cage traps (0.02 ± 0.02), while the Northern Brown Bandicoot had the highest (0.23 ± 0.03). Weekly detectability using the array of camera traps ranged from 0.02 ± 0.02 for the Delicate Mouse to 0.94 ± 0.06 for Short-eared Rock Wallaby. Of the 20 mammals, six were detected exclusively on camera traps (mainly macropods), two were not detected using this method, and 13 were detected by cameras and another method. Of these 13 mammals, five were detected with a higher probability on cameras deployed for one week than when using nighty live trapping methods.

Fire frequency was retained in 7% of bird, 3% of reptile, and 16% of mammal detectability models, with an exclusively negative effect on reptile and mammal detectability, and mixed effect on birds (Fig 3; S9–S11 Tables). Ruggedness was retained in 5% of bird, 10% of reptile, and 16% of mammal detectability models, with an exclusively negative effect on mammals, and a mixed effect on birds and reptiles. Time since fire was retained in relatively few detectability models.

Species richness

For mammals, our relative species richness map predicted higher richness in rugged terrain along the edge of the Arnhem plateau and in rugged areas in the south west, as well as across coastal lowlands in the north and west (Fig 5). This pattern largely reflects the influence of terrain ruggedness and vegetation cover on predicted occupancy for many of the mammals (Fig 3; S1 Fig). Relative reptile richness was predicted to be highest on the Arnhem plateau, and in fragmented areas of rugged terrain to the south west of the study region (Fig 5). In contrast, relative bird richness was highest across vast coastal and lowland areas due to the combination of several landscape and fire covariates (Figs 3 and 5; S1 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Relative species richness maps.

Richness maps for 20 mammals, 33 reptiles, and 84 birds across the Top End of northern Australia. Note, the scale bar is the summed probability of occupancy for each species, not the expected number of species, and richness refers to the pool of modelled species, not all species in the region.

https://doi.org/10.1371/journal.pone.0203304.g005

Influence of environmental covariates on occupancy

For mammals, fire and landscape variables were retained in a higher proportion of occupancy models than broad-scale climatic variables. It should be noted, however, that we fitted occupancy models to fewer mammals compared to other species groups, which meant that the relative contribution of covariates in percentage terms was highly sensitive to the choice of model for each mammal species. There is mounting evidence that contemporary mammal declines in northern Australia have been influenced by fire frequency [44] and extent [12]. Causal factors driving species distributions are not always evident from observed patterns of occupancy, and are more reliably demonstrated with studies of spatial dynamics (e.g., dynamic occupancy modelling [33]). Nonetheless, we observed the anticipated negative association between fire frequency and extent and occupancy for several mammal species [12,45]. However, positive associations also occurred for some mammals such as Antilopine Wallaroo and Black-footed Tree-rat. Higher occupancy in burnt areas may not necessarily be indicative of a species preference for those conditions, but an artefact of species being restricted to these regions due to other biotic processes not considered in our study. For example, Black-footed Tree-rat occur in lowland woodlands, a habitat type more prone to extensive fires, thereby contributing to the observed positive relationship between fire extent and occupancy.

Higher predicted mammal richness in rugged terrain followed by coastal lowlands to the north and west is in broad agreement with current hypotheses about where mammals are persisting in the Top End. Monitoring and distribution data indicate that range contractions have occurred from south to north, with earlier and more severe mammal declines inland [13], and greater resilience of populations of some species on Coburg Peninsula and offshore islands [46]. Previous studies concur with our findings of high mammal richness in rugged/rocky areas in the region [47]. Rocky upland areas contain a largely distinct mammal assemblage (e.g., rock rats and rock wallabies), with Arnhem Land Rock Rat and Black Wallaroo endemic to the Arnhem plateau. Mammals in rugged terrain may be able to persist longer as factors driving declines elsewhere may be diminished or absent in these areas. For example, rocky and rugged terrain presents barriers to fire, resulting in patchier and less damaging fires compared to less complex lowland areas where there are generally fewer impediments to fire spread [48]. While landscape variables were retained in many of our models, mammal declines are thought to be driven by more complex interactions between abiotic and biotic processes, such as altered fire regimes [12], predation by feral cats and/or grazing by feral herbivores [10,13], and the interaction of these factors [49,50], as well as interaction with rainfall and productivity [47]. Accurately accounting for and incorporating abiotic and biotic interactions into predictive models requires a good understanding of these ecological processes.

Unlike birds, most studies of reptile distributions across northern Australia have focused on the influence of habitat variables at finer scales than we used in our analysis [23,51], however, these studies do not attempt to predict distributions across landscapes. One study within Kakadu National Park considered 102 species and found that most were associated with a gradient of substrate and moisture availability, but rarely vegetation structure, with highest richness in rocky areas [23]. The results of our broader scale study concur with these findings by identifying high reptile richness in rugged areas–a useful proxy for rockiness–despite the fact that we only modelled a subset of the reptile community. Another study that sampled across the rainfall gradient of the Top End found that sites in high rainfall areas supported the most reptile species, and sites on clay soils supported fewer species [25]. Our modelling approach found little support for rainfall as a predictor of occupancy, however, clay content was retained in more reptile models than for mammals and birds, albeit still a small proportion.

In comparison to mammals and reptiles, relative bird richness is predicted to be highest across a large proportion of the study region due to the broad distributions of many birds throughout northern Australia [52]. The relatively high contribution of maximum temperature and fire frequency on bird occupancy is largely in agreement with a broad-scale study by Reside et al. [43], who found that temperature seasonality and rainfall was the greatest predictor of occurrence for 44 species across northern Australia, with a small, yet significant effect of fire frequency for some species. Several studies have found marginally higher diversity in higher rainfall areas across the region, by inferring distribution patterns of birds from site-based sampling [53]. Generalisations about the influence of rainfall and climate on species distributions are of limited value in isolation, because as we have demonstrated, occurrence may vary extensively within these gradients due to finer-scales factors such as soil, vegetation cover, and fire.

Detectability using traditional sampling methods

Unsurprisingly, we found high variability in the detection probability of birds, reptiles, and mammals using our standardised monitoring protocol. Contrary to our expectations, fire frequency had a negative influence on the detectability of mammals, reptiles, and birds, while ruggedness had both positive and negative effects on detectability. The role of site covariates such as these on detectability therefore requires further consideration using finer-scale covariate data. Camera trapping was relatively effective at detecting mammals compared with live trapping, which agrees with recent research for similar species [17]. Importantly, the effectiveness of camera trapping seemingly increased with mammal body size (e.g., Delicate Mouse 0.02 ± 0.02 versus Agile Wallaby 0.49 ± 0.06). This result highlights that camera traps should not be used exclusively if the goal is to detect ground-dwelling mammals of all sizes. However, given cameras can be left unattended at sites for long periods they may present a relatively cost-effective sampling method compared with live trapping, although cameras need to be retrieved and substantial effort is required to process photos. Evaluating the cost-effectiveness of competing approaches to surveying species, taking into account the trade-off between effort/cost and precision of the estimates of occupancy [54,55], would benefit the redesign of the Three Parks Program.

Limitations

We made several assumptions that may have influenced occupancy and detectability estimates. Firstly, the initial set of environmental covariates considered in the most complex models were identical for birds, reptiles and mammals. Alternatively, we could have specified unique models for each species group, or for each species, however, this would have increased the complexity of our analysis. Secondly, we ignored within year variability in detection probability or any differences in detectability between observers. Expanding our detectability models to test for observer and seasonal effects warrants further investigation. Thirdly, a mismatch in the initial resolution of covariate raster layers may have introduced uncertainty into the analysis. Finer-resolution mapping of the coarsely mapped covariates, particularly fire, will likely improve our understanding of their relationship with occupancy and detectability. Finally, the home range of species varied considerably relative to the area of surveyed sites, which influences our interpretation of occupancy. For species with home ranges similar in size to the sampling unit of 100 x 100 m (i.e. reptiles and small mammals), occupancy can be interpreted as the probability of a species being present at a site. However, many of the birds and large mammals will likely have home ranges much larger than this area, meaning occupancy for these species is better interpreted as probability of use. This assumes that the lack of closure between repeat surveys is random with respect to movement in and out of sites [33].

There are several caveats with our relative species richness maps. Firstly, species richness maps were generated by summing the mean predicted occupancy across species. This approach ignores the variance around each species occupancy estimate, assumes all species are given equal weight, and can potentially over-predict richness at species-poor sites and under-predict richness at species-rich sites [42]. Secondly, for numerous species we were limited to presenting only null (constant) models which did not contribute to richness maps, and models with single covariates or with associated low regression coefficients. These models may be deemed as having limited value for interpreting patterns of occupancy; however, they were retained as they present a starting point from which we can build upon as data becomes available. Thirdly, we could only fit models to 136 species, which is approximately half of the species recorded at sampling sites, and only a fraction of the ~600 terrestrial vertebrates known to occur in the region. This means our richness maps are biased towards the more common, widespread and detectable species, and do not represent the true richness of the study region. Very few threatened species contributed to our species richness maps because there were either too few detections or the models were deemed implausible (due to inadequate predictors or relatively few detections). Updating our models in the future using data from new or existing monitoring sites might improve the representation of threatened species in the species richness maps. Alternatively, occupancy maps for these species could be developed using alternative data (i.e. presence-only data) and modelling approaches. This is an interesting avenue of further research [56,57] but beyond the scope of this study.

We generated relative occupancy and species richness maps by predicting our models to environmental covariates mapped across the Top End. In doing so, we extrapolated our predictions outside the geographic area that was sampled. This is common practice in species distribution modelling, provided care is taken to avoid extrapolating to sites where the response of species to environmental covariates might differ [58]. We were cautious when extrapolating our predictions to not predict far outside the environmental domains (i.e., combination of predictor ranges) that were sampled. Furthermore, we did not predict occupancy to offshore islands as obvious biogeographical barriers may have prevented dispersal, and islands are geographically disjunct locations where species might respond differently to the environmental covariates studied. To assess how well our models matched the observed process, we ran goodness-of-fit tests (Sum of Squared Errors, Freeman-Tukey Chi-squared, and Pearson’s Chi-squared test) on the best models for all species. We could not use traditional measures of predictive performance, such as AUC, given our non-detections were not confirmed absences, due to imperfect detection [59]. Further research into goodness-of-fit testing and model validation is needed in the context of occupancy-detection models.

The utility of models as a planning tool in conservation is influenced by the selection of relevant predictors [39]. Our models were constrained by using only 11 predictors, and the absence of numerous others well known to drive species occurrence and distributions, but for which broad scale mapping is unavailable (e.g., cat occurrence, grazing pressure, habitat structure and complexity, etc.). In the absence of a definitive set of mapped covariates there are two options available to us when designing monitoring. (1) Not using a spatial model–this ignores not only the factors driving distribution that can’t be mapped, but all spatially mappable factors. Instead, we choose to use imperfect models to improve the design of monitoring acknowledging there are deficiencies. (2) Using mappable environmental predictors as ‘distal predictors’ [60] or surrogates of the actual proximal drivers of distribution. The use of ‘distal predictors’ requires prior knowledge of the interaction between biotic factors, or biotic and abiotic factors. For example, the abundance of feral cats or feral herbivores in some parts of northern Australia is known to be influenced by climate, topography and vegetation type [61]. Where these correlations exists, the estimated response of species to these environmental factors will incorporate these effects, albeit imperfectly, thereby acting as distal predictors of feral cats and stock in this example. When a map of cat and grazing pressure distribution is available at sufficiently high spatial resolution, it would make sense to incorporate them into the analysis.

Informing monitoring decisions

Our species-level occupancy and detectability modelling combined with our species richness maps will inform redesign of the Three Parks Program. Firstly, our occupancy and species richness maps can provide a starting point from which to assess how existing or proposed monitoring sites align with species distributions. Decision-makers can then consider how new sites could be positioned to maximise either the spatial coverage of species or the number of species detected at sites. Furthermore, the location of sites can be compared with the expected or known drivers of occupancy to ensure adequate sampling across gradients, and relocated if necessary to stratify more effectively, thereby improving the ability to study faunal responses those drivers. For example, comparison between covariate values at sampled sites with values across the entire study region in this study suggests that the existing sites are fairly representative of some covariates (S6 Fig; e.g., topographic ruggedness, clayiness, fire extent), but do not sample the full spectrum of others (e.g., fire patchiness, annual rainfall, vegetation cover)—warranting the need to truncate our predictions of occupancy and richness to the environmental domain for which we had data (see methods). Thus, a redesigned vertebrate monitoring program for the Top End might seek to improve monitoring in regions masked out of our occupancy predictions (in the south of the study region or on islands) to improve our understanding of species across these environmental gradients.

Our detectability estimates for the six sampling methods will help evaluate which methods are worth retaining in the revised Three Parks program. Detectability estimate will also inform how much sampling effort is required to detect species with a given level of confidence. The results presented here can also inform similar monitoring decisions for related species outside of the study region that lack rigorous detectability estimates. An added benefit of estimating occupancy and detectability is that the statistical power of a revised monitoring program can be assessed. Statistical power is related to the number of sites, initial occupancy, detectability, the effect size, and the Type I error rate [62]. Future work could use the occupancy and detectability values reported here to quantify how many sites, or how much survey effort is required to detect changes in species occupancy with sufficient power (e.g., [5,6]).

Our work demonstrates that high resolution, broad-scale occupancy-detection modelling can be achieved using repeat presence-absence data that accounts for imperfect detection. Modelling repeat presence-absence data is an advance on the more widely used presence-only models (commonly fitted with the Maxent software, e.g., [24,43]) for a range of reasons, including explicitly accounting for ‘false zero’ observations in survey data. This is important because ignoring false absences can lead to biased inferences about species distributions and over-confident precision around parameter estimates [33,63,64]. Importantly, accounting for imperfect detection also allows one to determine whether an observed difference in occupancy over space and time is because of changes in populations, or because of a change in the detection probability. Despite these benefits, only a small minority of studies measuring species richness, occupancy or range/distribution account for imperfect detection [2,4]. Our study demonstrates that this approach can be applied to a broad suite of species across relatively large spatial scales.