Co-occurring species often overlap in resource use and can interact in complex ways. However, shifts in environmental conditions or resource availability can lead to changes in patterns of species co-occurrence, which may be exacerbated by global escalation of human disturbances to ecosystems, including conservation directed alterations. We investigated the relative abundance and co-occurrence of two naturally sympatric mammal species following two forms of environmental disturbance: wildfire and introduced predator control. Using 14 years of abundance data from repeat surveys at long-term monitoring sites in south-eastern Australia, we examined the association between a marsupial, the common brushtail possum Trichosurus vulpecula, and a co-occurring native rodent, the bush rat Rattus fuscipes. We asked: Is the increase in abundance of common brushtail possums following control of an introduced predator associated with a decline in abundance of the bush rats?

Using Bayesian regression models, we tested hypotheses that the abundance of each species would vary with changes in environmental and disturbance variables, and that the negative association between bush rats and common brushtail possums was stronger than the association between bush rats and disturbance. Our analyses revealed that bush rat abundance varied greatly in relation to environmental and disturbance variables, whereas common brushtail possums showed relatively limited variation in response to the same variables. There was a negative association between common brushtail possums and bush rats, but this association was weaker than the initial decline and subsequent recovery of bush rats in response to wildfires.

Using co-occurrence analysis, we can readily infer negative relationships in abundance between co-occurring species, but to understand the impacts of such associations, and plan appropriate conservation measures, we require more information on interactions between the species and environmental variables. Co-occurrence can be a powerful and novel method to diagnose threats to communities and understand changes in ecosystem dynamics.

Study location

We used data from long-term annual monitoring that commenced in 2003 in Booderee National Park (BNP), Jervis Bay Territory, south-eastern Australia. BNP is on Indigenous land and is jointly managed by the Wreck Bay Aboriginal Community and Parks Australia. The 6600 ha park has a temperate climate, with an average annual rainfall of 1213 mm, spread evenly across the year (Bureau of Meteorology, 2021). Average temperatures range from 18.6-25.1°C in summer (January) to 9.9-16.1°C in winter (July) (Bureau of Meteorology, 2021). BNP supports a range of vegetation types such as heathlands, wetlands, forests, and woodlands. Two major fires have occurred in BNP over the past 20 years (in 2003 and 2017), with each burning approximately half of the park. A fox baiting program has been in place in BNP since 1999 and was intensified in 2003 to reduce the deleterious impact of this introduced predator on native prey species (Dexter et al., 2013).

Data collection

We surveyed small and medium-sized mammals annually each summer for 14 years at 109 permanent sites starting in 2003, with another 20 sites added in 2008. The sites were surveyed along 100 m transects with 2 large (30 x 30 x 60 cm) cage traps at the beginning and end of transects, small (20 x 20 x 50 cm) cage traps every 20 m between the large cage traps, and 10 Elliott traps (10 x 10 x 30 cm; Elliott Scientific Equipment, Upwey, Victoria) every 10 m (Figs 2) (Lindenmayer et al., 2008, 2016). Approximately 50% of the sites were surveyed each year (with the other 50% being surveyed the next year), depending on weather conditions (Lindenmayer et al., 2016). We recorded the number of individuals of both species caught at each site in a given year.

                We collected data on environmental and disturbance variables at each of our 129 sites. These data included visual estimates of the percentage of understorey and leaf litter cover in four 1 x 1 m subplots within 20 x 20 m survey plots during semi-annual vegetation surveys (Lindenmayer et al., 2008). We selected understorey and leaf litter as representative variables of the primary bush rat habitat, which are also components of common brushtail possum habitat (Callander, 2018; Cruz et al., 2012). We constructed a predictive model to fill the data gaps for those years when sites were not surveyed (MacGregor et al., 2020). We used monthly rainfall averages collected at the nearby Point Perpendicular weather station for the trapping period at each site (Bureau of Meteorology, 2021). We transformed both the vegetation variables and rainfall were transformed into quadratic functions using the poly() function in R (R Core Team, 2021).

                We used data on fire occurrence recorded on-ground since 2003, and fire history data collected by Booderee National Park over the past 50 years (Foster et al., 2017), specifically the number of years since the last fire at a site. To minimise possible inaccuracies stemming from incorrect fire dates, or the occurrence of unreported fires, we grouped the number of years since fire into 10-year blocks (i.e., 0-10, 11-20, 21-30, 30+ years).

Statistical analysis

We used Bayesian regression models with a hurdle step to test the response of species abundances to the selected variables using the brms package ver. 2.16.3 (Bürkner, 2017; 2108; Feng, 2021) implemented in R (R Core Team, 2021). These regression models employed Markov Chain Monte Carlo simulations, with four chains and a warm-up of 1000 iterations before sampling another 1000 iterations. We assessed model convergence by ensuring all Rhat values were <1.1 (Bürkner, 2017; 2108). The hurdle step consisted of two components: the first modelled the presence/absence of the response variable, and the second, conditional on the species being present, modelled the conditional abundance using a zero-truncated Poisson (Feng, 2021). We combined these two components in an analysis of unconditional abundance (Feng, 2021).

We created a regression model with bush rats as the response variable, and a regression model with common brushtail possums as the response variable. Both regression models included time, understorey cover, leaf litter cover and rainfall as covariates within the conditional abundance component. These variables were included to assess the variation in bush rat and common brushtail possum abundances with environmental variables (H1). Years since fire was included in the conditional abundance and hurdle step of both regression models as an explanatory variable, as it is a prominent disturbance within BNP, and past research has indicated that fire has a significant effect on small vertebrate populations (Arthur et al., 2012). The other species was also input into the conditional abundance and hurdle step of both regression models (i.e., common brushtail possums into the bush rat model, bush rats into the common brushtail possum model) as an explanatory variable to assess the co-occurrence effect between species (H2). We also included site as a random effect, and used the log of the number of Elliott traps as a control for the bush rat models, and the number of open cage traps as a control for the common brushtail possum models. The control variables account for varying trapping effort, and were selected based on the main trap-type that captures the relevant species (i.e., Elliott traps for bush rats, cage traps for common brushtail possums).

We performed a model selection procedure for both of the regression models, based on the selection for explanatory variables only. We chose not to perform model selection on the covariates (i.e., the environmental variables) as we were testing variation in species abundance in relation to the environment, and not predicting significant changes in abundance that we were with the co-occurrence and disturbance variables. Using model selection, we assessed the relevancy of our exploratory variables to changes in species abundance (H3). The chosen model was the most parsimonious, which was based on the simplest model which was within 2 leave-one-out cross validation (LOOIC) scores of the best fitting model.

We created ten variations of the regression models for each species, and assessed the fit of each variation using LOOIC (Tables 1 and 2) (Vehtari et al., 2017). LOOIC estimates the out-of-sample predictive fit by measuring the predictive accuracy for each data point using a variation of the expected log pointwise predictive density equation (Vehtari et al., 2017). LOOIC was selected as the appropriate method over other model selection methods as it is informative and was created for Bayesian models (Burnham & Anderson, 2002; Vehtari et al., 2017).