Occupancy and co-occurrence patterns of endemic mammals and introduced predators across a broad geographical gradient in eastern Australia

Invasive predators, land clearing and altered fire regimes have been implicated in species decline and extinction worldwide. Enhanced knowledge of how these factors interact and influence medium-sized mammals is warranted. We tested three hypotheses using occupancy data for a diverse mammal assemblage including three threatened species, five common species, two introduced mesopredators and an apex predator in eastern Australia. We hypothesised that occupancy of mammal species within the assemblage would be influenced by (i) the physical environment (rainfall, vegetation type and elevation), (ii) habitat disturbance (number of fires and habitat fragmentation) and (iii) mesopredator release, whereby occupancy and/or detection of medium-sized mammals are influenced by mesopredators, the feral cat (Felis catus) and the red fox (Vulpes vulpes), which are influenced by an apex predator, the dingo (Canis familiaris). We utilised camera-trapping data from 173 sites (692 camera locations) across a north–south gradient spanning ~ 1500 km in eastern Australia. Although hypotheses i (physical environment) and ii (habitat disturbance) are not mutually exclusive, we show that the variables considered in each were only weakly correlated. We conducted occupancy modelling to investigate the physical environment and habitat disturbance hypotheses. We conducted co-occurrence modelling to investigate interactions between species. The physical environment hypothesis best supported occupancy models for six mammal species: red-necked pademelon (Thylogale thetis), bandicoots (Isoodon macrourus and Perameles nasuta), swamp wallaby (Wallabia bicolor), red-necked wallaby (Macropus rufogriseus), eastern grey kangaroo (Macropus giganteus) and feral cat. The disturbance hypothesis best supported occupancy models for four mammal species: long-nosed potoroo (Potorous tridactylus), red-necked pademelon and both mesopredators. Support for the mesopredator release hypothesis was equivocal. Large macropods showed site avoidance towards the red fox. Four species showed higher detection at sites where mesopredators were not detected. The fox showed a negative detection interaction to the dingo and the cat did not. Our study highlights how factors such as rainfall, land clearing, elevation and number of fires influence the occupancy of species within a diverse mammal assemblage at the macroecological scale. Our findings have implications for the conservation of threatened species in managed landscapes and suggestions for further research following the recent 2019–2020 wildfires.


Introduction
The world is facing a biodiversity crisis with the continuing extinction and decline of species globally Rosenberg et al. 2019). While human activities have increased the rate of species extinctions to more than 1000 times the natural background rate (Pimm et al. 1995), managed conservation areas are experiencing secondary legacy effects from human activities long after initial threats have subsided. Threats such as habitat disturbance, introduced predators and altered fire regimes impact biodiversity within protected areas and these factors are likely to be exacerbated outside of conservation lands. The plight of many ground dwelling mammals is such that their persistence requires rewilding efforts where populations are translocated to "safe-havens" such as fenced enclosures and islands where introduced predators and land clearing are absent . However, gaining a better understanding of how threats and environmental factors influence species persistence outside of such safe havens will be ongoing and important for biodiversity conservation.
Australia's diverse assemblage of terrestrial mammals is amongst the most unique and distinctive in the world (Holt et al. 2013). In the past two centuries Australia's mammals have experienced a higher rate of extinction compared to the mammal fauna of any other nation (Woinarski et al. 2015). Thirty-four of Australia's 273 endemic mammal species have become extinct in the past 200 years (Woinarski et al. 2019) and at present, there are 30 species listed as endangered under Federal legislation with a further 46 listed as vulnerable (Environment Protection and Biodiversity Conservation Act 1999). Among those particularly affected are medium-sized terrestrial mammals (Burbidge and McKenzie 1989;Cardillo and Bromham 2001). Although the broad factors that drive global extinctions are well known (Gibbons et al. 2000;Stuart et al. 2004;Skerrat et al. 2007;Szabo et al. 2012;Sánchez-Bayo and Wyckhuys 2019), further investigation of factors causing declines in managed landscapes are paramount to inform conservation and management and arrest further declines.
The single largest threat to Australian mammals at the continental scale is predation by introduced mesopredators, primarily the European red fox (Vulpes vulpes) and feral cat (Felis catus) (Short and Smith 1994;Kinnear et al. 2002;Woinarski et al. 2015;Radford et al. 2018). At smaller spatial scales, major threats implicated in declines include habitat loss and fragmentation (Law and Dickman 1998;McAlpine and Eyre 2002;McAlpine et al. 2006;Reside et al. 2019) and altered fire regimes (Woinarski et al. 2015). Recent studies suggest that there may be a synergistic effect between the presence of introduced predators and disturbance factors such as fire (Robley et al. 2016;Hradsk et al. 2017). It is unknown how introduced mesopredators respond to disturbance factors associated with landscape fragmentation, however, it is suspected that they may benefit (Catling and Burt 1995a).
Apex predators are widely acknowledged to shape ecosystem processes through direct lethal effects that influence prey populations (Ripple et al. 2014) and through non-lethal effects that alter prey species behaviour (Shrader et al. 2008;Laundré et al. 2014). Apex predators may also influence species assemblages indirectly through the lethal and nonlethal effects they have on mesopredators (Gordon et al. 2015). Understanding the role of apex predators is fundamental to maintaining biodiversity in terrestrial ecosystems (Ritchie and Johnson 2009). However, maintaining or restoring the ecological function of apex predators is often contentious due to the direct threat they pose to humans and threatened species (Allen and Fleming 2012;Augusteyn et al. 2020) and due to predation on livestock and associated economic impacts (Fleming et al. 2006;Muhly and Musiani 2009;Chetri et al. 2019). Indeed, these conflicts with human interests frequently lead to population control of apex predators (van Eeden et al. 2018).
The mesic forests of south-eastern Australia provide refugia for many mammalian taxa that are congeneric with species that have become extinct or suffered major range contractions in arid ecosystems. The challenges of understanding the factors that influence mammal populations over this large geographic region is the presence of overlapping gradients in environmental factors that may influence populations and ultimately species persistence. For example, dingoes (Canis familiaris) and mesopredators are widespread and subject to varying levels of population control but there are pronounced gradients in topography, rainfall and temperature that may also influence endemic mammal populations. In addition, disturbances which influence mammal populations such as wildfire and land-clearing are widespread Bradstock 2010;Reside et al. 2019). Consequently, largescale studies that span these environmental gradients are needed to compliment fine-scale approaches as they encompass much wider gradients and heterogeneity in such factors.
In this study, we use data from "Wildcount" a government funded large geographicscale camera-trapping program. We investigate three broad hypotheses concerning the relative influence that 'bottom-up' factors (the physical environment and habitat disturbance) and 'top-down' factors (predators) have on occupancy patterns of medium-large (1.3 to 19 kg) mammals ( Fig. 1), including three threatened species, in diverse habitats of eastern New South Wales (Tables 1 and 2). Our first hypothesis was that many of our study species will be influenced by factors arising from the physical environment, namely, rainfall, elevation and/or vegetation type (Catling and Burt 1995b;Claridge and Barry 2000;McHugh et al. 2019). We predict that rainfall will drive site productivity and this may enable higher reproduction and therefore persistence of species. In contrast, elevation may also represent site productivity but it will also be associated with topographic complexity which may facilitate persistence (e.g. McDonald et al. 2017). Vegetation type was also investigated as it reflects species food and shelter requirements.
Our second hypothesis, which is not mutually exclusive of the physical environment hypothesis, was that habitat disturbances such as fire history and habitat fragmentation/ land clearing will drive patterns in the occupancy of mammals (Claridge and Barry 2000;Norton et al. 2015). The disturbance factors we chose to investigate were the number of fires over 50 years, and several metrics relating to habitat fragmentation; perimeter/area ratio (P/A ratio) of the reserve surrounding each study site, and the percentage of cleared land surrounding each site with buffers of 5 km for endemic mammals and 10 km for introduced predators.
Our third hypothesis was that mesopredators would have a strong influence on mammal assemblages and that occupancy of mesopredators in turn would be influenced by the presence of dingoes, the apex predator in our system (Colman et al. 2014;Hunter et al. 2018;Rees et al. 2019). We tested our hypotheses in two steps. First, we conducted multiseason occupancy modelling to investigate whether any site-based covariates relating to the physical environment or habitat disturbance can explain variation in the detection of our 1 3 mammal species. Second, we conducted co-occurrence modelling to investigate interactions between each of the predators (dingo, fox, cat) and threatened and non-threatened mammals. We also investigated whether the three predators showed patterns of co-occurrence that are dependent or independent of each other.

Study area
This study was conducted across the eastern portion of New South Wales (NSW) in south-eastern Australia (Fig. 2). Marked climatic differences across the study area reflect changes in latitude and elevation. The coast extends 1460 km north to south. Elevation ranged from 0 m on the coast up to 1809 m on the tablelands west of the coast. The highest rainfall occurred in sub-tropical northern NSW (annual average = 2912 mm) and generally decreased along a latitudinal gradient south-ward, with lower rainfall in the temperate south-east corner (annual average = 1523 mm). The lowest rainfall occurred in the southwestern slopes region (annual average 509-675 mm).
The study area contained a range of different vegetation types with site numbers proportional to their extent: dry sclerophyll forest was present at 102 sites, wet sclerophyll forest was present at 29 sites, grassy woodland was present at five sites), heathland communities was present at seven sites, semi-arid was present at seven sites, rainforest communities were present at four sites, wetlands were present at six sites and alpine complex were present at two sites.

Survey design
Wildcount is a state-wide fauna monitoring program established by the New South Wales Office of Environment and Heritage in 2012 (EPA 2018). The project covers a large spatial gradient across eastern NSW (Fig. 2) that includes a broad range of habitats, climatic zones and fire histories. Approximately 200 monitoring sites have been established across 146 conservation reserves. Sites were selected using a stratified random selection of 1 × 1-km grid cells overlayed by a 20 × 20-km grid to prevent site clustering and to maximise dispersion. Sites were considered candidates for inclusion if they (i) occurred on NSW reserves and (ii) were intersected by a 4WD access track. Each site contained four camera traps (Reconyx HyperFireTM PC 800 (Infra-red flash) spaced within a 1 × 1-km grid approximately 500 m apart, which were surveyed annually over a 2-week period. Each camera was fixed to a tree at a height of approximately 1 m, directed at a stainless-steel tea strainer which contained a mixture of peanut butter and oats as a lure. The lure was positioned at a distance of 2 m from the camera trap at a height of approximately 20 cm. Cameras operated for 24 h per day and were set to take 3 images per trigger with no quiet time between triggers.
Images were identified to species by trained staff and trained volunteers using an image tagging system in the program Exifpro. Images were assigned an identification confidence   Table 2 Testing of hypotheses  of interaction between predators  and potential prey detection The column headings show the direction of difference for predictions if the predators have an influence on the prey. The parameter estimates are shown for each species pair and whether the values are the same (=) or different (> or <). The cells in bold highlight estimates consistent with the detection hypothesis p A = Probability of detecting species A (predator) given only species A is present, r A = probability of detecting species A (predator) when both species A and species B (subordinate species) are present, r B/A = probability of detecting species B (subordinate species) given that both species are present and species A (predator) is also detected, r B/a = probability of detecting species B (subordinate species) given that both species are present and species A (predator) is not detected, p B = probability of detecting species B (subordinate species) given that only species B is present

Landscape covariates
We selected six landscape covariates hypothesised to influence the occupancy of our target species. These were vegetation type, fire history (number of fires), rainfall, elevation, perimeter/area ratio of conservation reserves and % foliage projective cover (FPC) of woody vegetation within a 5-km and 10-km radius surrounding sites. We developed

Bioregions
Bioregions are large geographic areas characterised by common geophysical, hydrological and climatological features that influence the distribution of plant communities, faunal assemblages and ecosystem processes at large spatial scales (Thackway and Cresswell 1995). Under the interim biogeographic Regionalisation of Australia (iBRA) framework, 85 bioregions are recognised of which 17 occur across NSW. The Wildcount dataset cover nine of these Bioregions (Fig. 3A). Many of our species do not occur uniformly across the broader study area. Attempting to model species where they do not occur is likely to produce spurious results by inflating the sample size of non-detection sites. Therefore, we restricted our species data to the Bioregions where species predominantly occurred.

Vegetation
We

Rainfall
Rainfall is highly variable between years in eastern Australia due to the El Niño-Southern Oscillation (ENSO) and it is well appreciated that rainfall can drive substantial fluctuations in the abundance of forest mammals. For this reason, we used rainfall in the 3 years preceding each year of camera trapping. We retrieved gridded (5 × 5-km) average monthly rainfall data across Australia from the Australian Bureau of Meteorology (www. bom. gov. au/ jsp/ awap/ rain/ archi ve) in March 2018. Using ArcGIS, monthly rainfall data were firstly merged into three years of average rainfall and then each year was merged to form an average rainfall raster for the survey years. The raster was then clipped to the state of NSW and continuous values of average rainfall were extracted to each monitoring site within the study area (Fig. 3C).

Elevation
We retrieved digital elevation data from Geoscience Australia (www. ga. gov. au/ scien tifictopics/ natio nal-locat ion-infor mation/ digit al-eleva tion-data). The Digital Elevation Model (DEM) were derived from data generated from the 1 s Shuttle Radar Topography Mission (SRTM) acquired by NASA in February 2000. We extracted continuous elevational data to each camera monitoring site within the study area using ArcGIS (Fig. 3D). We derived an elevation value for each site by calculating the average of elevation from the four camera trap locations.

Fire history
We retrieved fire history data from OEH (http:// datas ets. seed. nsw. gov. au/ datas et) in August 2016 (Fig. 3E). The data comprise a number of attributes that relate to each fire boundary polygon including fire type (prescribed burn or wildfire), year of fire, area and perimeter of fire. Due to high correlation between these fire variables, we retrieved a single fire history covariate: number of fires that have occurred over a since 1970 which ranged between zero and nine.

Land clearing/habitat fragmentation metrics
We retrieved woody vegetation and Foliage Projective Cover (FPC) spatial data from the NSW Department of Planning and Industry (DPIE) (https:// datas ets. seed. nsw. gov. au/ datas et/ lands at-woody-extent-and-folia ge-proje ctive-cover-fpc-ver-2-1-25m-20087 355d). These data are derived from a time series of Landsat FPC images from 1988 to 2008 and use time series decision tree statistics to identify areas of woody vegetation and assign FPC values to them at a resolution of 25 m. We generated circular buffers around the NSW Wildcount sites at 5 km and 10 km from which we extracted FPC values. FPC values were then divided into two discrete categories (i) non-woody vegetation and (ii) woody vegetation. These values were used to generate a percentage value for cleared (grassy) land within each buffer area as a surrogate for cleared/farm land (Fig. 3F). We also generated a perimeter area ratio covariate based on the perimeter and area of each conservation reserve that Wildcount sites were located within. This was achieved by calculating the perimeter and area of reserves containing survey sites and dividing the perimeter by the area.

Single species occupancy modelling
We conducted multi-season occupancy modelling to test the physical environment and disturbance hypotheses. Modelling was implemented in program PRESENCE version 9.3 (USGS Patuxent Wildlife Research Centre, Laurel MD, USA). Multi-season models include four parameters representing different probabilities: initial occupancy (psi), detection (p), colonisation (gamma) and local extinction (epsilon) (MacKenzie et al. 2003). Our primary focus was on covariates that influenced species occupancy so we allowed colonisation and local extinction to be constant over time. For three species, gamma converged on zero so this parameter was fixed to zero. We allowed detection to either remain constant over time or to be year specific. We conducted occupancy modelling for 12 species (Table 1; Fig. 1) including several recognised as threatened: the long-nosed potoroo, red-legged pademelon and Parma wallaby. We constructed weekly detection histories showing whether each species was detected (1) or not (0) across the three 2-week monitoring periods.
We firstly compared models for the probability of detection with occupancy set to constant. We then retained the top detection model and examined the influence of covariates, representing our hypotheses, on the probability of occupancy. Covariates used were drawn from the physical environment and the disturbance variables described above. Models were ranked using Akaike's Information Criterion (AIC) from lowest to highest. The difference in AIC (∆AIC) was calculated between the top model and other models. A model was considered most plausible in explaining the data if it differed to the next model by > 2.0 ∆AIC (Burnham and Anderson 2004). Increasing values of ∆AIC > 4.0 indicated models had poor support (Burnham and Anderson 2004).
We commenced by fitting univariate models that aligned with our hypotheses and proceeded in the following sequence of covariates: (i) elevation, (ii) rainfall, (iii) number of fires, (iv) vegetation type, (v) cleared vegetation within a 5-km buffers for endemic mammals and 10-km buffers for predators surrounding sites and (vi) perimeter to area ratio for conservation reserves surrounding sites. We constructed dummy variables (0 or 1) to represent each vegetation type comprising the four forest and woodland types and uncommon types pooled into an 'other' group. This was conducted due to the way Presence reads variables within the design matrix. Rainfall, elevation and number of fires were standardised by subtracting their mean and dividing by their standard deviation. We then fitted models with two or more covariates by including those from models in which model weight (w) was > 0.1 (see Eyre 2007). We did not assess model fit because methods to assess this have not yet been devised for multi-season models (MacKenzie et al. 2018).

Co-occurrence modelling
We employed the mesopredator release hypothesis to help explain co-occurrence patterns within our occupancy data. Co-occurrence modelling can be used to investigate the interactive effects between pairs of species (MacKenzie et al. 2018). We constructed stacked detection histories of paired species, including each of the predator species with subordinate species, and with each other. We used Presence version 9.3 (USGS Patuxent Wildlife Research Centre) to model the species pairs using the formulation of Richmond et al. (2010). In this formulation one species is designated as the dominant species (A) and the other as a subordinate species (B). We designated a predator as the dominant species and the potential prey as the subordinate species. The probability of occupancy is estimated for the dominant species (psiA) and the subordinate species when the dominant species is present (psiB/A) and absent (psiB/a). The probability of detection is estimated for each species when the other is absent (pA, pB), when each of the species is present and detected (rA, rB/A), and for the subordinate species, when the dominant species is present and not detected (rB/a).
Co-occurrence modelling followed a two-step process. We firstly fitted models to investigate influences on occupancy and then retained the top model and investigated 1 3 the influence of changing the detection parameters. We used landscape covariates for individual species identified as influential from single-species modelling. We compared models with and without the covariates previously identified to be influential. Covariates were retained only where they improved model fit by > 2.0AIC (see below). We then investigated whether the prey species occurred at sites independently of the predator (psiB/A ≠ psiB/a). To investigate the detection parameters all detection parameters were initially estimated separately (pA, pB, rA, rB/A, rB/a). We then investigated whether detection of the predator differed when the prey species was present or not present (rA ≠ pA); whether the detection of the prey at sites with the predator differed depending on whether the predator was detected or not (rB/A ≠ rB/a); and whether detection of the prey differed when the predator was present or not (pB ≠ rB/A). Whether the estimated parameters that were compared were equal or different was used to provide an initial finding about interactions between species pairs. The modelling produces two derived parameters, phi (the species interaction factor) and delta (the detection interaction factor) (MacKenzie et al. 2018). Phi represents interaction between species as a function of occupancy, whilst delta estimates interaction between species as a function of detection. These parameters do not demonstrate that interaction is occurring but that occupancy or detection of two species is not random with respect to each other. If phi or delta are > 1.0, it suggests a positive association between the two species (targeting of prey by a predator or that both prefer the same habitat); if phi or delta are < 1.0, it suggests a negative association (avoidance). If phi or delta are ~ 1.0 then species associate independently. If the 95% confidence intervals of these parameters overlapped 1.0 we viewed this as lack of clear evidence of an interaction. The estimation of delta in Presence does not account for the difference between pB and rB so we also rely on whether pB, rB/A and rB/a differ to infer a detection response to the predators.
Co-occurrence models were ranked using AIC. The top model was treated as the most plausible if the second model differed by > 2.0AIC. Where the second model was within 2.0AIC of the top model we still used the top model to estimate parameters because it had fewer parameters and the choice made little difference to the parameter estimates. If individual models did not converge they were removed.

Long-nosed potoroo
The long-nosed potoroo was detected on 46 occasions across 13 sites and had a naïve occupancy of 0.12. Detection was estimated as constant across occasions (0.84 ± 0.07). The most plausible occupancy model was one that included the number of fires (Table 3). Based on model weight this model had 2.9 times more support than the null model. Gamma was estimated at 0.03 ± 0.01 and epsilon at 0.22 ± 0.12. The probability of potoroo occupancy decreased with an increasing number of fires (Fig. 4B).

Red-legged pademelon
The red-legged pademelon was detected on 66 occasions across 11 sites and had a naïve occupancy of 0.14. Detection of the red-legged pademelon was estimated as constant across occasions (0.77 ± 0.08). The most plausible occupancy model was the null model (Table 3). The addition of covariates showed no reduction in AIC. Gamma was estimated at 0 and epsilon at 0.26 ± 0.11.

Red-necked pademelon
The red-necked pademelon was detected on 80 occasions across 20 sites and had a naïve occupancy of 0.26. Detection was estimated as constant across occasions (0.81 ± 0.04).
The most plausible occupancy model was one that included cleared vegetation within a 5-km radius, rainfall and number of fires (Table 3). This model had 3.5 times more support than the next model. Gamma was estimated at 0 and epsilon at 0.25 ± 0.08. The probability of occupancy decreased with increasing number of fires (Fig. 4B), increased with increasing rainfall (Fig. 4A) and increased with increasing percentage of cleared land (Fig. 4B).

Parma wallaby
The Parma wallaby was detected on 29 occasions across 11 sites and had a naïve occupancy of 0.14. Detection was estimated as constant across survey occasions (0.68 ± 0.10). The top occupancy model included the elevation covariate and was   (Table 3). In the null model initial occupancy was estimated as 0.09 ± 0.04. Gamma was estimated at 0.03 ± 0.02 and epsilon at 0.42 ± 0.14.

Bandicoots
Bandicoots were detected on 333 occasions across 103 sites and had a naïve occupancy of 0.59. Long-nosed bandicoots and northern brown bandicoots were pooled due to uncertainty in distinguishing the two species in camera trapping imagery. Detection was estimated as constant across survey occasions (0.63 ± 0.03). The most plausible occupancy model included annual rainfall (Table 3). This model had a model weight of 1.0. Gamma was estimated at 0.10 ± 0.03 and epsilon at 0.15 ± 0.04. The probability of occupancy increased with increasing rainfall (Fig. 4A).

Red-necked wallaby
The red-necked wallaby was detected on 253 occasions across 70 sites and had a naïve occupancy of 0.41. Detection was estimated as constant across occasions (0.75 ± 0.03). The top occupancy model included rainfall and elevation (Table 3). This model had a model weight of 0.99. Gamma was estimated at 0.06 ± 0.02 and epsilon at 0.17 ± 0.04. The probability of occupancy decreased with increasing rainfall but increased as elevation increased (Fig. 4A).

Swamp wallaby
The  (Table 3). This model had a model weight of 1.0. Gamma was estimated at 0.34 ± 0.09 and epsilon at 0.10 ± 0.03.

Eastern grey kangaroo
The eastern grey kangaroo was detected on 246 occasions across 72 sites and had a naïve occupancy of 0.42. Detection was estimated as constant across occasions (0.63 ± 0.03).
The most plausible occupancy model included rainfall (Table 3). This model had a model weight of 1.0. Gamma was estimated at 0.02 ± 0.02 and epsilon at 0.09 ± 0.03. The probability of eastern grey kangaroo occupancy decreased with increasing rainfall (Fig. 4A).

Feral cat
The feral cat was detected on 105 occasions across 66 sites and had a naïve occupancy of 0.38. Detection was estimated as constant across occasions (0.19 ± 0.02). The top model contained elevation and Perimeter/Area ratio (Table 3). This model had a model weight of 0.98. Gamma and epsilon were fixed at 0.0. The probability of occupancy increased as elevation increased (Fig. 4A) and decreased as the Perimeter/Area ratio of conservation reserves increased (Fig. 4B).

Red fox
The red fox was detected on 287 occasions across 95 sites and had a naïve occupancy of 0.55. Detectability was estimated as constant across occasions (0.63 ± 0.03). The top occupancy model included the percentage of cleared land within a 10-km radius (Table 3). This model had a model weight of 1.0. Gamma was estimated at 0.16 ± 0.03 and epsilon at 0.12 ± 0.04. The probability of occupancy increased as the percentage of cleared land increased (Fig. 4B).

3 Dingo
The dingo was detected on 34 occasions across 25 sites and had a naïve occupancy of 0.14. Detectability was estimated as constant across occasions (0.30 ± 0.11). The top occupancy model included woodland vegetation however this model was equally plausible when compared to the null model (Table 3). Gamma was estimated at 0.08 ± 0.03 and epsilon at 0.74 ± 0.16.

Summary of hypothesis evaluation
The physical environment hypothesis had support for most generalist species whereas the habitat disturbance hypothesis was supported for sensitive habitat specialists and introduced mesopredators. Of the covariates representing the physical environment hypothesis, rainfall had the strongest overall influence over endemic mammal occupancy with two small forest dwelling species (red-necked pademelon and bandicoots) showing a positive response to rainfall and three herbivorous large macropod species (eastern grey kangaroo, red-necked wallaby and swamp wallaby) showing a negative response. Vegetation type was not influential over any species and elevation was influential over the red-necked wallaby and the feral cat. Of the covariates representing the habitat disturbance hypothesis, the long-nosed potoroo and red-necked pademelon were negatively influenced by number of fires and the red-necked pademelon was negatively influenced by % of cleared land in 5-km radius of sites. The two introduced mesopredators benefited from habitat disturbance where the red fox showed a very strong response to % cleared land in 10-km radius of sites and the feral cat showed a strong negative response to the perimeter area ratio of reserves.

Red fox with other species
All the top models except that for the potoroo, the red-necked wallaby and feral cat included separate parameters for psi B/A and psi B/a (Table 4). However, only the swamp wallaby and the eastern grey kangaroo appeared to show site avoidance of the red fox as indicated by values of phi < 1.0 (Fig. 5). The red-legged pademelon, Parma wallaby and bandicoots showed positive site associations with the fox however only the red-legged pademelon showed phi values > 1.0 (Fig. 5). Four species including the red-legged pademelon, Parma wallaby, red-necked pademelon and bandicoots showed much higher detections at sites without the fox (p B ) compared to sites with the fox (r B/A ) ( Table 4). Estimates of delta for these species overlapped or were approximately 1.0 suggesting no detection interactions where both species occurred (Fig. 6).

Feral cat with other species
The parameters for psi B/A and psi B/a were estimated as equal in the top models for potoroos and foxes (Table 5). None of the estimates of phi suggested avoidance of cats with values at or above 1.0 (Fig. 5). The red-legged pademelon and eastern grey kangaroo showed phi Table 4 Comparison of 2-species models for the red fox (species A) and subordinate species (species B) Only the top two models or those within 2AIC of the top model ψ A = probability of occupancy for species A, ψ BA = probability of detecting species B given species A is present, ψ Ba = probability of detecting species B given species A is absent, p A = Probability of detecting species A given only species A is present, r A = probability of detecting species A given both species A and species B are present, r B/A = probability of detecting species B given that both species are present and species A is also detected, r B/a = probability of detecting species B given that both species are present and species A is not detected, p B = probability of detecting species B given that only species B is present. E = elevation; R = rainfall, F = number of fires, 10-km = % of cleared vegetation in 10-km buffer around site, 5-km = % of cleared vegetation in 5-km buffer around site values above 1.0 suggesting positive a site interaction. None of the species showed a detection avoidance at sites where cats occurred though red-legged pademelon, Parma wallaby, bandicoots and red-necked wallaby had higher detection at sites where cats did not occur ( Table 5). The estimates of delta for the red-legged pademelon, and bandicoots showed a positive association with cat detection (Fig. 6).

Dingo with other species
The parameters for psi B/A and psi B/a were estimated as equal in the top models for longnosed potoroo, red-legged pademelon, swamp wallaby, eastern grey kangaroo and the red fox suggesting there was no site interaction with these species (Table 6). Three species (Parma wallaby, red-necked pademelon and red-necked wallaby) showed higher occupancy at sites occupied by dingoes compared to where dingoes were absent with Parma wallaby and red-necked wallaby showing phi values > 1 suggesting positive site interactions   Table 6 Comparison of 2-species models for the dingo (species A) and subordinate species (species B) Only the top two models or those within 2AIC of the top model are shown ψ A = probability of occupancy for species A, ψ BA = probability of detecting species B given species A is present, ψ Ba = probability of detecting species B given species A is absent, p A = Probability of detecting species A given only species A is present, r A = probability of detecting species A given both species A and species B are present, r B/A = probability of detecting species B given that both species are present and species A is also detected, r B/a = probability of detecting species B given that both species are present and species A is not detected, p B = probability of detecting species B given that only species B is present. R = rainfall, F = number of fires, 10-km = % of cleared vegetation in 10-km buffer around site, 5-km = % of cleared vegetation in 5-km buffer around site  5). Bandicoots showed a negative site association with the dingo (Fig. 5). Several species showed detection responses to the dingo (Table 6). The long-nosed potoroo and eastern grey kangaroo had higher detection when dingoes weren't detected whereas Parma wallaby and red-necked wallaby had higher detection at sites where dingoes did not occur. At sites where foxes and dingoes occurred foxes showed lower detection when dingoes were detected. Both the eastern grey kangaroo and red fox had delta estimates < 1.0 suggesting avoidance in the presence of dingoes (Fig. 6).

Discussion
We investigated the extent to which three non-exclusive hypotheses can explain occupancy patterns of medium and large endemic mammals and their introduced predators across a broad geographic scale in eastern Australia. We did not expect uniform responses by species to the factors representing each hypothesis due to differences in species' life histories. Our results supported the physical environment hypothesis for six mammal species (Table 7), the habitat disturbance hypothesis was supported for four mammal species (Table 7) and support for the mesopredator release hypothesis was equivocal (Table 7). Findings from our large-scale study present both opportunities and constraints which can guide further research that aims to untangle relationships between our endemic mammal species, introduced predators and the biotic and abiotic environment. The varied responses of species to our hypotheses outlines some important findings regarding single-species occupancy modelling over the broad geographic gradient in this study and also outlines  some challenges with conducting co-occurrence occupancy modelling over such gradients. We now discuss each hypothesis in turn and conclude with implications for further research, conservation and management.

Physical environment hypothesis
Our study area encompassed a very broad geographic scale so we expected that environmental factors would have a strong influence on occupancy of most of our study species. Thus, it was not surprising that rainfall was an influential predictor of mammal occupancy. Rainfall is a well-known factor driving geographic patterns in mammal assemblages (Heaney 2001;Olff et al. 2002). Forest specialists such as bandicoots and red-necked pademelons showed a positive response to rainfall while grassland and open habitat specialists such as eastern grey kangaroos, and red-necked wallabies were negatively influenced by rainfall. We equate this response to a broad influence on plant community type that reflects the rainfall gradient. The lack of influence of rainfall on the three threatened species suggests that productivity may not be a primary determinant of their persistence in the landscape.
Although elevation can have a pronounced influence on species' distributions (e.g. Ramírez-Bautista and Williams 2019; Campera et al. 2020) we found only one endemic species, the red-necked wallaby, responded (positively) to a broad gradient (1800 m) in elevation. The lack of response from other endemic species may reflect that our sampling was mostly focused on forest habitats and that most of our focal species are habitat generalists. The Parma wallaby generally occupies areas of high elevation across our study area, however we only found a weak effect of elevation on this species which may have been hampered by a small sample size. One important relationship we identified was that the probability of occupancy of the feral cat increased with elevation. The mechanism behind this relationship may be that it is a reflection of higher productivity at high-elevation sites that induces higher abundance of prey such as small mammals (Bateman et al. 2010). Small mammals comprise approximately 70% of feral cat diet in Australia (Murphy et al. 2019) and high elevation sites where small mammals are abundant may support higher number of cats.

Habitat disturbance hypothesis
We hypothesised that some mammal species would be influenced by fires, percentage of cleared land and the perimeter/area ratio of reserves. Fire is a frequent disturbance in Australia's forested ecosystems and plays a significant role in altering the composition and structure of vegetation communities and habitat (Bradstock 2010). We found a negative relationship between the number of fires over the past 50 years and the occupancy of the threatened long-nosed potoroo and the non-threatened, red-necked pademelon. This is a novel finding for the red-necked pademelon but not for the long-nosed potoroo. For the long-nosed potoroo, the model that included number of fires showed a modest improvement to model fit compared to the null model. However, previous studies that have investigated the post-fire response of the long-nosed potoroo suggest that occupancy will increase over time following wildfire, with > 15-20 years being optimal, as this fire interval allows for the development of dense ground cover habitat (Claridge and Barry 2000;Norton et al. 2015). The red-necked pademelon showed a stronger negative response to number of wildfires when compared to the long-nosed potoroo which may have resulted from a larger sample size. The ability of species to recover post-fire when vegetation structure is redeveloping will depend on species mobility and the availability of unburnt refugia in the surrounding landscape (e.g. Leonard et al. 2014). Recent studies in temperate Australia suggest that there may be a synergistic effect between fires and introduced predators, whereby fire implemented in the presence of introduced predators have negative effects on mediumsized mammals including potoroos (Robley et al. 2016;Hradsk et al. 2017). However, this finding may not be uniform across the distribution of this species. Long-nosed potoroo and predator activity did not change following prescribed burns in sub-tropical Australia (McHugh et al. 2020). Greater caution may be required in using fire in this species' habitat in the southern but not northern part of its distribution due to the varied abundance of foxes.
One limitation of our study was that we were not able to ascertain the response of species to prescribed burns due to a smaller sample of sites where prescribed burns had occurred. Wildfires will generally encompass large areas and have a high intensity/severity (Leonard et al. 2014;Chia et al. 2016), whereas prescribed burns will generally be small in scale and conducted during cool weather with low-moderate intensity/severity (McHugh et al. 2020). Consequently, wildfires can be expected to be more detrimental to potoroos compared to prescribed burns. Of particular concern is that the 2019-2020 wildfires in eastern Australia overlapped approximately 40% of the potoroo's NSW location records (DPIE 2020).
Habitat loss and fragmentation due to clearing of vegetation is considered a threat to biodiversity at a global scale (Pimm and Raven 2000). Native habitats in eastern Australia have been extensively cleared historically ) and this continues in contemporary times (Reside et al. 2019). Understanding the influence of this on biodiversity is critical for conservation. We identified a negative response by the red-necked pademelon to the extent of clearing within a 5-km radius of survey site. Other native species were not influenced perhaps because our survey sites were contained within conservation reserves. Some important findings relating to habitat disturbance are that we found that introduced predators were influenced by land clearing/habitat fragmentation covariates; red fox occupancy increased as clearing increased within a 10-km buffer of sites, whilst feral cat occupancy decreased as the P/A ratio of conservation reserves increased. These are important findings to help understand factors that influence these introduced predators in fragmented landscapes of eastern Australia. These findings can guide conservation efforts targeting these species regarding population control.

Mesopredator release hypothesis
We predicted that dingoes would have a negative association with foxes and cats, which in turn would have a positive association with their prey. We found no evidence of a site-scale interaction between dingoes and feral cats or red foxes. However, we found a negative detection interaction whereby fox detection was significantly lower when dingoes were also detected compared to when they were absent. This suggests that foxes show temporal avoidance of dingoes. Previous studies have indicated that smaller predators may avoid larger predators due to the risk of intra-guild predation (Palomares and Caro 1999;Linnell and Strand 2000;Brawata and Neeman 2011). Lack of site interaction between dingoes and feral cats is consistent with studies finding little support for the idea that cats modify their use of habitat in response to dingo activity (McHugh et al. 2019;Fancourt et al. 2019).

3
There is an emerging body of literature that is questioning the inference that can be drawn from multi-species occupancy modelling with regard to survey design i.e. spatiotemporal extent, types of detectors used, covariates used, choice of field methods and statistical tools (Devarajen et al. 2020). Experimental evidence and observation support the notion that species interactions influence co-occurrence patterns, however Blanchet et al. (2020) question to what extent the signal of interaction can be retrieved from observational data. A limitation in the present study is that we used weekly detection periods to improve the detection of co-occurrence signals which in-turn may weaken the inference that can be drawn, a challenge that presents itself where dominant and subordinate subject species have limited spatial overlap (Letnic et al. 2009;Brawata and Neeman 2011;Colman et al. 2014Colman et al. , 2015. For example, the detection window of one week may provide opportunity for species to co-occur without directly experiencing each other through using risk averse behaviours. Finding a trade-off between opening the detection window and deriving strong inference from field data remains a challenge for co-occurrence modelling and we hope that through identifying this challenge in the present study, further research may benefit when considering survey design and statistical methods (Devarajen et al. 2020).
In addition, covariates that may be influential but are not considered may weaken inference drawn from co-occurrence modelling (Devarajen et al. 2020). For example, habitat preference may mediate co-occurrence of species (Estevo et al. 2017;McHugh et al. 2019) and many species in this study show a preference towards habitat structure at finer scales (McHugh et al. 2019). Fine scale habitat structure data were not available across the broad scale of this study though it is likely to influence several species in this study and potentially mediate co-occurrence patterns. In support of our findings that did not find a site interaction between the dingo and the red fox, previous studies investigating dingo-fox interactions have found that fox abundance/activity tends to be lower at sites where dingo populations are not controlled and that there is little spatial overlap between them (Letnic et al. 2009;Brawata and Neeman 2011;Colman et al. 2014Colman et al. , 2015. We found no evidence that foxes avoid dingoes at the site scale, however, this result may be because there was relatively little spatial overlap between these species in our dataset. The naïve occupancy of the dingo was 0.14 whilst that of the red fox was 0.55. The relatively low occupancy of dingoes may reflect ongoing control of populations over much of the study region (Colman et al. 2014;Ballard et al. 2020). Consistent with the idea that there is limited spatial overlap between dingoes and foxes, a previous study which deployed cameras at 298 sites across nine forested conservation reserves in NSW detected dingoes at 18% of sites but detected foxes at just 2% of sites (McHugh et al. 2019).
We found no evidence of a negative site association of the mesopredators and the smaller macropods and bandicoots. Foxes had a positive site-based association with the red-legged pademelon and a negative association with larger macropods (swamp wallaby and eastern grey kangaroo). An explanation for this may be that small macropods have the ability to co-exist with foxes through displaying risk averse behaviours i.e. concealment within dense habitat (Signorell et al. 2010) whereas larger macropods cannot and are more inclined to use rapid movement to evade predation (Creel et al. 2005). Furthermore, detection of bandicoots and small macropods was much higher at sites without foxes compared to sites with foxes (Table 2), supporting the idea that behavioural avoidance may be taking place.
Similarly, the feral cat had positive site-based occupancy associations with several species (red-legged pademelon, Parma wallaby and red-necked wallaby) and all of these species had higher detection values at sites without cats (Table 2). Dingoes showed a positive site-based occupancy association with some small macropods (Parma wallaby and red-necked pademelon) and a negative association with one large macropod: the eastern grey kangaroo, however, bandicoots also displayed a negative response to sites where dingoes were present (Fig. 5). Detection avoidance of the dingo was displayed by the eastern grey kangaroo which is frequently preyed upon by the dingo (Wright 1993) and interestingly this species showed a positive detection interaction with the smaller red fox. Previous literature suggests that red foxes pose a significant threat to eastern grey kangaroo joeys that can limit population growth of kangaroos (Banks et al. 2000) though its likely adult eastern grey kangaroos can avoid fox predation. Despite the apparent co-existence of endemic mammals and introduced mesopredators within the mesic conservation reserves of eastern Australia, concern arises of the impact of feral cats in high-elevation areas (Fig. 4A) and foxes in fragmented areas (Fig. 4B) where their occupancy is relatively high and therefore the threat to endemic species is elevated.

Implications for further research, conservation and management
Our study has provided knowledge and insights on the factors the drive occupancy for endemic mammal species and introduced predators. Rainfall is a key driver of endemic mammal occupancy in eastern Australia. The wetter regions of the study area support habitat for small macropods including threatened species; Parma wallaby, red-legged pademelon, long-nosed potoroo. Although these species have distributions that surpass the NSW state border, high rainfall regions eastern Australia will be important in providing habitat refugium for these threatened species and also can guide decision making with regard to prioritising recovery sites.
Habitat disturbance remains an insidious threat to threatened mammal species in eastern Australia, in particular the threat of wildfires. During late 2019 and early 2020 south-east Australia experienced some of the largest wildfires on record. The 2019-2020 Australian bushfire season occurred during an extended period of low rainfall and record-breaking temperatures (Filkov et al. 2020). It is estimated that the mega-fires burnt ~ 97,000 km 2 across southern and eastern Australia, which is considered habitat for 832 species of native vertebrate fauna (Ward et al. 2020). The present study demonstrated the negative influence of number of fires on the long-nosed potoroo and the red-necked pademelon. The 2019-2020 wildfires in eastern Australia overlapped approximately 40% of the potoroo's NSW location records (DPIE 2020) which is alarming given the impacts that introduced predators may have in temperate regions on long-nosed potoroo distribution following fire (Robley et al. 2016;Hradsk et al. 2017), and we recommend that immediate and further research is required on long-nosed potoroo populations effected by recent wildfires.
In addition, the present study found that landscape fragmentation directly influences introduced mesopredators: the red fox and feral cat across the study area. Our finding that occupancy of the red fox increased with percentage of cleared land is consistent with a study conducted 25 years previously (Catling and Burt 1995a). Given the paucity of fox records in forested areas of north-east NSW and the diversity and abundance of mammal species (Catling and Burt 1997), this region can clearly be identified as refugium for many threatened mammals in addition to the long-nosed potoroo, Parma wallaby and the red-legged pademelon. Threatened species in fragmented south-east temperate regions will benefit from ongoing efforts to control fox populations through lethal baiting.
There is a view that mesic ecosystems of eastern Australia support lower feral cat densities compared to their xeric counterparts ) and therefore the impact of cats on biodiversity may be lower in mesic ecosystems (Doherty et al. 2017;Woinarski et al. 2017). However, our finding that feral cat occupancy increased with elevation has not been identified previously and raises concern that cats may potentially have more severe impacts on fauna in forested high-elevation areas than previously recognised. We suspect that the mechanism that allows for such high cat occupancy in high-elevation forested areas across our study region is that they provide ample food resources for feral cats (small mammals, reptiles and birds) and also habitat structural complexity that allows cats to avoid potential predators. Given that invasive mesopredators have had devastating impacts on biodiversity in xeric systems (Burbidge and McKenzie 1989;Cardillo and Bromham 2001;Woinarski et al. 2015) we recommend that further research is conducted on the ecology and impacts of cats in high elevation forests.
Our study and its design was not without its limitations, including that lure type used was suboptimal for detecting predators (Hanke and Dickman 2013) and our camera traps were located off road within vegetated habitat rather than on tracks, which may result in suboptimal detection of predators (Read et al. 2015;Geyle et al. 2020). However, lure type was standardised across all sites and each site had four widely spaced camera traps, which undoubtedly increased detectability of predators (Stokeld et al. 2015). Also, sites were intersected by a 4WD access track and therefore roads/tracks were in close proximity to camera trap locations. Further modelling on predator interactions within this mammal assemblage may benefit from considering these limitations.
The present study has outlined some important findings regarding single-species occupancy modelling over a broad geographic gradient in eastern Australia and also some challenges associated with co-occurrence occupancy modelling. The very broad study area gives rise to a wide range of variables that could influence our mammal assemblage and by necessity, some of our predictor variables are coarse. Also, alternative hypotheses may be influential at finer scales when considering single-species occupancy modelling, for example habitat structure and resource availability. Similarly, alternative hypotheses may be considered for two-species co-occurrence interactions (Blanchet et al. 2020). However, the benefit of this large-scale approach is that it has evaluated mammal responses to predictor variables over a range of values that is often not possible in studies conducted over smaller geographic scales, and it has encompassed many species not possible at finer scales. Further considerations regarding co-occurrence occupancy modelling should be given to the spatiotemporal overlap of species which will remain a challenge when considering the dingo and the red fox.