Estimating biocontrol agent spread: A case study using introduced dung beetles

Insect biological control agents may be released in relatively few locations; thus, the rate of spread is critical to predicting the spatial distribution of the agent. Exotic dung beetles have been introduced into Australia since the 1960s to bury dung, reduce bush fly populations, increase pasture productivity and improve livestock health. We use a stochastic cellular automata model to predict the historical occupancy and abundance of introduced dung beetles in the South West agricultural region of Western Australia. The model predicts the spread of four species of dung beetles in six‐month time steps for 10 years following initial release. Our model includes release locations, species‐specific ecological parameters, dung resource availability and weather variables. The average rate of spread varied between species from a high of ca. 79 km/year (Euoniticellus intermedius) to a low of 28 km/year (Onthophagus taurus). Thus, after 10 years, the area of spread varied between 99,175 and 34,175 km2 and abundance varied from 88,100 to 4,189 beetles/km2. These findings provide an estimate of the spread patterns of dung beetles. The model can be used to guide future dung beetle release programmes in Australia and elsewhere.


INTRODUCTION
Biocontrol agents are usually introduced to 'control' a weed or pest species.That is, reduce population size to reduce economic losses (Barratt et al., 2018;Caltagirone & Doutt, 2003;Conti et al., 2021).Two characteristics of a successful biocontrol agent are rapid spread across its target range and effective control of the pest or weed.For instance, a planthopper introduced as biocontrol agent against smooth cordgrass in Washington, U.S., spread 200 m over 3.5 months (Grevstad et al., 2003).Understanding the spread of biocontrol agents is crucial to assess the likelihood of establishment, and the ecological and economic impacts of a biocontrol programme (Heimpel & Asplen, 2011).
There are several reasons for predicting the rate of spread of biocontrol agents.First, modifying management practices to allow survival of agents (Barratt et al., 2018), for example, reducing pesticides (Bielza et al., 2020).Second, assessing agent efficacy postestablishment.Third, for the maintenance of the agent population, as the agent population may change, even die off and require reintroduction after extirpation (Grevstad, 1999).Thus, understanding spread can assist in designing environmental management strategies, guiding economic impact assessments, supporting introduction decisions (McNitt et al., 2019;Stansbury et al., 2002;Tonini et al., 2013Tonini et al., , 2014) ) and improving programme delivery.Determining spread rate allows control targets to be set, and conversely to support population growth of the agent.Despite these advantages, our understanding of spread is dominated largely by weed biocontrol agents (insects and fungi) (Culliney, 2005;Paynter & Bellgard, 2011;Thomas & Reid, 2007).
Many factors may influence rate of spread: climate, ecological interactions (Bebber et al., 2014), propagule size (Duncan, 2016;Grevstad et al., 2012) and resources.Most spread models focus on climate from the native range; however, the transferability of climate envelope models for introduced species is highly debated (Duncan et al., 2009;Liu et al., 2020), as they do not include biotic interactions and resource availability (Duncan et al., 2009).Ecological interactions such the presence of predators, competitors and diseases can also change the rate of spread (Crowder et al., 2019).Therefore, it is important to include factors that may change spread patterns for a more accurate spread model of introduced species.
Spread models are useful when monitoring is absent.Monitoring biocontrol agents is essential for the assessment of ecological and economic impacts, however it is rarely done (Schaffner et al., 2020) due to expense.Biocontrol programmes are expensive (often decades long), and usually, funds are not available for post-release.Field monitoring is expensive (Wainwright & Mulligan, 2013), especially in remote areas.So modelling may be useful when monitoring programmes are not feasible, as models can provide rapid results and reliable predictions when using high-quality data (Drew et al., 2011;Wainwright & Mulligan, 2013).
Proof of the usefulness of spread models comes from their use for accidental pest introductions.These include dengue fever mosquito (Kumar Acharya et al., 2018), Colorado beetle (Pulatov et al., 2014), bronze bug (Montemayor et al., 2015) and spotted wing drosophila (Ørsted & Ørsted, 2019).Different model approaches include two-dimensional kernels to predict pine beetle spread (Koch et al., 2020;Strohm et al., 2016), and a random forest model using satellite images of tree cover, temperature and rainfall data to predict spruce beetle spread (Hart et al., 2017).Such models have had practical applications in ecology, conservation biology (see Rossiter-Rachor et al., 2023), and pest management (e.g., Karnal bunt and mountain pine beetle) and are an important tool in management programmes to assess ecological impacts and economic impacts (Celesti-Grapow & Ricotta, 2021;Cooke & Carroll, 2017;Stansbury et al., 2002).
Exotic dung beetles have been used as biological control agents in many countries, including Australia, North America and New Zealand to bury livestock dung (Duncan, 2016;Forgie et al., 2018;Pokhrel et al., 2021).By burying dung, dung beetles prevent dung accumulation and control fly infestations in livestock systems (Edwards, 2007a;Hughes, 1975;Pokhrel et al., 2021), recycle nutrients, improve plant growth and reduce nematode parasites affecting ruminant livestock (Badenhorst et al., 2018;Bang et al., 2005;Doube, 2018;Johnson et al., 2016).They also provide benefits indirectly to society by reducing pollution related to dung contamination of pasture and water runoff (Brown et al., 2010) and supressing nuisance flies.The value of dung beetles to the cattle industry has been estimated to be up to £367 million each year in the United Kingdom (Beynon et al., 2015).In Mexico, the estimated annual value of dung beetle services to the cattle industry was up to US $423.6 per animal unit (Lopez-Collado et al., 2017).
A total of 52 dung beetle species were imported into Australia between 1969and 1984(Tyndale-Biscoe, 1996).There was five years of monitoring, mostly to determine if the beetles had established viable populations (Ridsdill-Smith & Hall, 1984).After that, there were no further monitoring attempts to predict the spread of dung beetles across the wider landscape.The lack of extensive monitoring, especially post-initial releases, hinders our understanding of impacts of this biocontrol programme and spread patterns.
In this article, we modelled the spread of dung beetles introduced to the South West agricultural region of Western Australia (WA).Our model is a single-species model applied to four dung beetle species and estimated species occupancy and abundance in six-month time steps.We applied the model to four species known to establish in WA: Onthophagus binodis [released 1972], Euoniticellus intermedius [1972], Onitis alexis [1973] and Onthophagus taurus [1975] (Edwards, 2007b;Matthiessen et al., 2014;Ridsdill-Smith & Hall, 1984).We aimed to estimate the rate of spread, the abundance and the area occupied by each species 10 years after release.We chose 10 years because it encompassed the time period with detailed release records and was three times the average time of establishment (3.5 years) for exotic dung beetle species in Australia (Edwards, 2007a).

Dung beetle release data
We retrieved release data containing information on species name, release locations, quantity released, month and year of the releases from Tyndale-Biscoe (1996) (Figure 1.) and used them as the starting points of spread in our model (see Table S1 footnote for initial release details).We converted the month and year of the releases into halfyear intervals starting in the first half (January-June) of 1972, when the first dung beetle species, Ont.binodis, was initially released (see Table S1 footnote).We converted release locations from geographic (latitude and longitude) into projected coordinates (EPSG:3577).

Stochastic cellular automata model
We predicted dung beetle occupancy and abundance using a stochastic cellular automata model.These models simulate spatial units as grid cells and can represent dynamic processes (Pitt et al., 2009), such as the spread of wildfires (Jiang et al., 2020;Josephson et al., 2020) and the spread of invasive species across a landscape (Tonini et al., 2014).Stochastic cellular automata can represent large areas due to their computational efficiency (Tonini et al., 2014).We based our model on the termite spread model developed by Tonini et al. (2014).We extended this model by incorporating resource availability, dung beetle species-specific biological parameters and spatially varying adult survival probability.We ran the model using parallel processing, which is a method to execute tasks simultaneously on multiple processors to improve model speed.A diagram of the data utilised for this model is given in Figure 2.
We used a Monte Carlo (MC) technique with multiple simulation runs to account for the uncertainty in the predictions, and results are expressed as the mean with 95% CI of all the simulations.The R code for the spread model is provided in Appendix S3.Model parameters for each species are shown in Table 1.
The model simulated the spatial-temporal changes in the occupancy and abundance of dung beetles.We used 6-month time step, as dung beetle life or reproductive cycles are usually 6 or 12 months long.For each time step, beetles in cells increase in age, and those exceeding the 'MaxAge' parameter die.After each cell has been updated, a number of new offspring are generated based on the number of reproductive females in the cell and the 'MaxNOffspring' and 'EggToAdultSurvival' parameters.The dispersal of offspring to new cells is then simulated with a two-parameter dispersal kernel (Clark, 1998): where Γ() is the gamma function, α is the average flight distance parameter ('FlightDistance'), c is the shape parameter and x is the distance from parent cell.In the present study, c was set to 1, making the distribution exponentially bounded.
Once offspring were generated, the resulting number of reproductive females in each cell was calculated by randomly sampling from a lookup table of pre-computed simulations (for computational efficiency).These simulations calculated the number of reproductive females by assigning offspring a random position within a cell.
Assuming a sex ratio of 0.5 (Martínez et al., 2019), female offspring within a mating attraction distance of 1,000 m from a male were assumed to be reproductive based on the study by Beaudoin-Ollivier et al. (2003), thus increasing the abundance of reproductive females in that grid cell by 1.We assumed that the maximum number of females a male can find and successfully mate was 30.We constrained beetle numbers in two ways.We capped the maximum number of beetles supported by resources in a cell using a carrying capacity parameter (CCapacity) and randomly remove adult beetles depending on the 'ProbAdultSurvival' and 'MaxSurvival' parameters.We describe the calculation of these parameters in the sections below.

Carrying capacity
We defined carrying capacity as the maximum number of individuals that the environment can support, considering habitat availability, amount of dung available and dung beetle food requirements.In our model, the carrying capacity estimate limits the number of reproductive females per grid cell.For model simplicity, we treat carrying capacity as a fixed parameter for each grid cell (no changes in the livestock number or pasture area).

Habitat
Habitat availability was given by the pasture area, as it contains dung.
We retrieved the earliest available digital land use (based on satellite imagery) data from ABARES (2022) for 1992-1993.We defined habitat in grids as 'grazing natural vegetation', 'grazing modified pastures' and 'irrigated modified pastures' as suitable for dung beetles, with a grid cell value of 1.All other land use categories were defined as nonsuitable and set to zero (Figure 1).We increased the land use cell size by approximately 1.11 km Â 1.11 km (0.01 ) (ABARES, 2022) to 5 km Â 5 km grid cell size using the nearest neighbour interpolation method.

Total dung production
We define food availability as the total amount of dung produced Local Government Area (LGA) (Figure S1).We chose 1972 as the base year because this was when the first 'successful' dung beetle species (Ont.binodis) was introduced to the intensive agriculture area of Western Australia.The ABS data reported the number of livestock for cattle and sheep per LGA.
Dung availability is measured by: where D s denotes the amount of dung produced (kg day  et al., 1995;DataGene, 2002;Lewer et al., 1992, and weight from: Browning et al., 1995;Handcock et al., 2019;Lowman et al., 1993;Primary Industries Standing Committee, 2007).To match the spatial resolution of our model, the total dung production (Figure S1) per LGA was distributed spatially to grid cells (Figure S2) through dasymetric mapping using the census 1981 shapefile (Australian Bureau of Statistics, 2018).Dasymetric mapping is used to convert vector (polygon) data to raster (gridded) data.

Species dung requirement
The amount of dung in kg required for survival, development and reproduction of an individual of each dung beetle species, i, is defined by the following: where R i is the amount of dung in kg required for the survival, development and reproduction of an individual.This is defined by the species fecundity, F i (i.e., number of brood-balls produced by an actively reproductive female per time step), and multiplied by the species brood-ball weight (B i Þ in kg (see Table S3).A brood-ball is defined as the buried mass of dung containing one egg (Halffter & Edmonds, 1982) or multiple eggs (Klemperer, 1982).Adult beetles feed on the liquid dung components to survive and develop and use the solid dung to reproduce and build brood-balls (Ridsdill-Smith et al., 1982).For this reason, we used fresh brood-ball weight B i ð Þ for our calculations.We assumed that the liquid component of the brood-ball is enough for survival and development of an individual, whereas the solid component will be used for reproduction.
In addition to the carrying capacity per grid cell, we also present the average number of individuals that the available dung in kilograms can support, by dividing it by the species brood-ball weight We estimated the carrying capacity per LGA in each time step: Posteriorly, the total carrying capacity per Local Government Area (C s,i Þ was distributed spatially to the suitable areas (Figure 3) using the dasymetric mapping technique.Equation ( 5) estimates the carrying capacity (maximum number of individuals) per grid cell: The carrying capacity per grid cell is defined by G s and is estimated by dividing the total carrying capacity per LGA (C s,i Þ by the total number of grid cells that are suitable habitats for dung beetles, H.The suitable habitat includes grid cells that are categorised as 'grazing natural vegetation', 'grazing modified pastures' and 'irrigated modified pastures'.

Adult survival probability
To estimate the adult survival probability by grid cell, we employed the maximum entropy MaxEnt algorithm (version 3.1.;Phillips et al., 2004Phillips et al., , 2006)).MaxEnt models use sample points and environmental layers as model inputs to provide probability scores.The probability scores describe how likely a species will be found at a location, within or without of its native range (Phillips et al., 2006).
We retrieved species presence points from the Global Biodiversity Information Facility (GBIF) portal (https://www.gbif.org/)for the four dung beetle species up to 2020.To define the species native range, we cleaned the data (code provided in Appendix S3) using information on species origins from Edwards (2007b) and removed outliers using the 'CoordinateCleaner' R package (Zizka et al., 2019).
Sampling bias is a common issue in species distribution modelling, due to uneven sampling effort (Phillips et al., 2006).To reduce sampling bias, we kept only one occurrence record within each grid cell.Spatial filtering has been found to reduce overfitting and improve model performance (Boria et al., 2014) et al., 2011).
For model optimisation, we used the 'ENMevaluate' function of the "ENMeval" R package (Muscarella et al., 2014) to assist in the choice of MaxEnt 'default' settings.We converted these raster outputs to 5 kmÂ5 km grid cells containing values ranging from zero to one, where zero means no survival probability and one means optimal survival probability.The adult survival probability was then used as an input in the model to assist in the prediction of each species occupancy and abundance.We removed individuals that fell into a grid cell where the adult survival probability is close to zero.
We interpreted the outputs of MaxEnt as 'ProbAdultSurvival'; the MaxEnt probabilities were then used to assign an adult survival probability to each grid cell for each species.At any one time point, the survival of adults is the product of MaxEnt (ProbAdultSurvival) and the maximum survival probability of an adult reproductive female (MaxSurvival).

Model development through sensitivity analysis and parameter calibration
Our overall approach to model development is presented in Figure 3.
The first step was to conduct sensitivity analyses to assess the importance of model parameters (Pannell, 1997;Wainwright & Mulligan, 2013).We did this to identify the model parameter that would likely have the most influence on the model results, which could then be tuned using model calibration to improve model performance (see below).To define the parameter with the most influence on the model results, we conducted a sensitivity analysis by varying the following parameters: 'EggToAdultSurvival', 'FlightDistance', 'Max-NAlates' and 'ProbAdultSurvival'.As the structure of the model was the same for all species, we conducted the sensitivity analysis for E. intermedius only.We ran 10 MC simulation for 10 time steps by varying the four parameters one at a time at three levels (i.e., À20%, 'default', +20%) while keeping the other parameter values fixed.The 'default' values are the 'best-bet' parameter values, and the selected levels are within a reasonable range.For the 'ProbAdultSurvival', every single 'default' grid cell was multiplied by À20% and +20% using the 'raster calculator' tool in QGIS (QGIS.org,2022).
After determining the most influential parameter, we used a model calibration procedure to find the value of this parameter that provided the best fit of the model output to the field data.The field data used for model calibration were monitoring data collected in the study region from 1977 to 1981 (see Figure S3A for locations and further details of the data used).Model calibration is the process of testing the model outputs by comparing these with field data, and then adjusting the parameters to improve the fit between predictions and observations (Mulligan & Wainwright, 2013).The choice of the parameter that will be used for model calibration depends on the results of the sensitivity analysis (Mulligan & Wainwright, 2013): the parameter most sensitive to changes is the target of model calibration.
Our model calibration procedure involved running the spread model for the period 1977-1981 (matched to observed data) and varying the value of the most sensitive parameter from À50% to +50% in steps of 10% while holding all other parameters constant.
We did this for three species but not Ont.taurus as it was not present in the observed data.We ran these models for 100 MC simulations.
This number of MC simulations provided sufficient precision at a low computational cost, similar to Tonini et al. (2013).After running each model, we used the 'zonal statistics' tool in QGIS (QGIS.org,2022) to calculate the mean predicted occupancy values for each six-month time step for each of the dung beetle monitoring locations (each with a 10 km buffer to account for coordinate inaccuracies).The model output for a raster cell ranges from 0 (0% of MC simulation runs) to 1 (100% of MC simulation runs), which we interpreted as probabilities of species presence in a raster cell at a time step.These probabilities were transformed to class predictions by assigning a '1' (i.e. a predicted presence) to average raster values above 0.5 and a '0' (i.e. a predicted absence) to average raster values less than or equal to 0.5.We then assessed our model predictions (presence/absence at a site during a time step) against the field data.A correct classification rate was used as a measure of fit to help us determine the best-fitted model.The parameter value for the most sensitive parameter was selected as the one that gave the highest correct classification rate score, and this value was used in the final model.

Model validation
We ran the spread model with 1000 MC simulations per time step for all species using the final parameter values selected by model calibration.
We assessed the performance of the models by comparing model predictions to monitoring data collected in the study region from 1990 to 1994 (Figure S3B, Dadour & Cook 1997).As these data extended beyond our analysis period, we ran an extended model of 45 years to validate our model.To do so, we ran 10 MC simulations per time step for 45 years (90 time steps) for the three species present in the 1990-1994 data.We tested the model predictions and calculated the correct classification rate for each model using the model calibration procedures.The 1990-1994 field dataset did not contain records for one of our species, E. intermedius.
We therefore also assessed the model predictions against the field data collected from 1977 to 1981.F I G U R E 3 Process for the sensitivity analysis, parameter calibration and validation for the model using data from (a) Ridsdill-Smith and Hall ( 1984) and (b) Dadour and Cook (1997).

Sensitivity analysis and parameter calibration
The sensitivity analysis identified the most sensitive parameter as 'FlightDistance', with variation in this parameter resulting in the greatest variation in the area occupied (ca.50 to 270,000 km 2 (Figure 4).
The value of the 'FlightDistance' was calibrated for species with sufficient data: Ont.binodis, E. intermedius and Oni.alexis (refer to Figure 3 and Table S6).A 'FlightDistance' value of minus 10% of the default value from Fincher et al., 1983 (i.e., 4 km per year or 2000 m per time step), predicted the highest correct classification rate for all the three species: correctly predicting 63.16% of the spatial-temporal spread for Ont.binodis, 78.94% for E. intermedius and 100% for Oni.alexis.
Thus, the selected 'FlightDistance' value was 1800 m per time step (2000 m Â 0.9).See Table S6 for more detail of the calibration results.

Model validation
Model validation showed around 80% correct classification across both periods (Table 2).Overall, false positives are more common in our models than false negatives.The model validation for 1977-1981 had the highest performance for Oni.alexis, while the model validation for the period of 1991-1994 had the highest performance for Ont.taurus.

Stochastic cellular automata model
The estimated values for the 'ProbAdultSurvival' and 'CCapacity' parameters are given in Table S7 Our findings indicated a rapid rate of spread of E. intermedius of approximately 79 km/year, which is consistent with reports from a few years after its initial release in Australia (Seymour, 1980) In contrast, our findings diverge from the existing literature on the spread of Ont.taurus.Our model predicted a spread rate of ca.28 km/year in WA, while in the United States, a rate of 275 km/ year was reported (Fincher et al., 1983).However, Ont. taurus was accidentally introduced in the United States, so the first record in 1971 (Fincher et al., 1983) may not be accurate.This accidental introduction was likely from just one location, whereas the deliberate introduction in WA were from many locations, likely influencing spread rates in WA.
The modelled spread results for dung beetles in WA are within the range of previous beetle spread models and field studies using re-capture data.A spread model predicted a spread rate of 79 km/year for the mountain pine beetle in Alberta, Canada (Cooke & Carroll, 2017).Field studies estimated the rate of spread for Digitonthophagus gazella to range from 34 km/year to 808 km/year from the United States to Mexico (Barbero & Lopez-Guerrero, 1992;Kohlmann, 1994).This wide variation in spread rates may be due to re-capture studies emphasising the farthest distances travelled.The mean dispersal distances are often not reported, even though it is more useful for spread models (Grevstad et al., 2003).This high variation in spread rates is reported for various dung beetle species (see   Hall, 1984) and model validation using field data from 1990 to 1994 (Dadour & Cook, 1997).These limited monitoring data constrain the testing of model predictions comprehensively; however, our model results do show satisfactory predictions of introduced dung beetles (70.3%-97.3%correct predictions at 45 years).Model parameter calibration and validation with observed field data are fundamental to improve model accuracy and performance (Wainwright & Mulligan, 2013).Our model predicts the species spread for a period of 10 years there is model uncertainty relating to long-term spread predictions, especially because of the uncertainty in climatic changes such as extreme events over time (Cooke & Carroll, 2017) and changes in habitat structure (Powell & Bentz, 2014).Even though  2020) highlight the benefits of spread models as a complementary tool to ground monitoring surveys.However, field data collection is not always possible due to cost, logistics and time constraints (Wainwright & Mulligan, 2013).Spread models are a useful analysis tool for planning and monitoring biocontrol programmes.

Appendix S1 and
They can provide detailed information to assist stakeholders in evaluating the time to meet the aims of species introductions.
Our model can assist stakeholders to understand the time before a desired level of dung burial or fly reduction is achieved.Data on the occupancy and abundance of introduced species are necessary to assess programme effectiveness, implement management strategies and allocate financial resources (Cheney et al., 2018).Such models can also facilitate programme economic valuation (i.e., economic benefits related to the agent introduction) and serve as an additional tool for monitoring programmes when the budget for sampling is limited and in remote areas where field sample monitoring is challenging.

CONCLUSIONS
This article investigated the rate of spread and abundance of introduced dung beetle species using resource availability, climatic variables and species-specific ecological parameters as model inputs.
Spread models are commonly used to predict the rate of spread of pest species.However, in our article, we demonstrate the importance of predicting spread patterns for decisions related to biocontrol programmes.Our model can be used by government agencies as a tool for predicting the rate of spread of introduced dung beetle, and to guide future introduction by targeting locations where species are not found.In addition, the model can be used by landholders as a low-cost alternative to monitor biocontrol agent post-release, in cases where ground-truthing data are limited or budget is limited for field surveys.
Future work is required for a more detailed assessment of the spatialtemporal spread by incorporating species competition.This research will also benefit from an economic evaluation to better understand the historical impacts of this biocontrol programme on the Australian livestock sector.

F
I G U R E 1 Release locations of Ont.binodis (purple square), E. intermedius (black circle), Oni.alexis (green diamond) and Ont.taurus (red triangle) in the South West agricultural region (black) of WA (white) between 1969 and 1984 sourced from Tyndale-Biscoe (1996).Only points containing information on the number of individuals released and release date are presented in this figure.The boundary for the South West agricultural region was retrieved from the Department of Primary Industries and Regional Development (2018).Habitat suitability grid cells (5 kmÂ5 km) are defined as suitable if the land use categories are 'grazing modified pastures' (light grey), 'grazing natural vegetation' (dark grey) or 'irrigated modified pastures' (blue); otherwise they are classified as non-suitable (white).A black line indicates the boundary of the South West agricultural region of WA.
by cattle and sheep.We use 1971-1972 data retrieved from the Australian Bureau of Statistics (ABS) (Australian Bureau of Statistics, 1972.As the base year for the number of animals per We defined C s,i as the carrying capacity by LGA calculated from dung produced (kg day À1 ) by LGA (D s Þ divided by the amount of dung required for survival, development and reproduction of an individual of each dung beetle species (R i Þ and then multiplied by the number of days in one time step 365 2 À Á . Strohm et al. (2016) developed a partial differential equation spread model for mountain pine beetles and found areas with high density of healthy trees-food resources for the beetle-heavily determined the spread dynamics.Similarly, our model results showed a greater beetle density in areas with a high density of livestock density and dung.Future research could include temporal changes in resource availability to improve our understanding of how changes in stocking density may affect dung beetle spread.Monitoring programmes can assist modellers in identifying the model constraints and serve as a guide to improve model performance and accuracy.However, there was only limited monitoring of dung beetles in the study area in the few decades post-release, thus limited empirical data of their spatial distribution.We conducted model calibrations using field data from 1977 to 1981 (Ridsdill-Smith & T A B L E 2 Model validation results derived for the model outputs (occupancy) and compared with Ridsdill-Smith and Hall (1984) (n = 57) and Dadour and Cook (1997) (n = 111) presence/absence records.

F I U R E 5
Predicted abundance (number of individuals per grid cell) in the South West agricultural region of Western Australia (Department of Primary Industries and Regional Development, 2018) after 5 years (left) and 10 years (right) from each species initial release.there are uncertainties, model predictions can be improved by including the main drivers involved in the spread dynamics of each species (Cooke & Carroll, 2017) such as competition.A recommendation for future work is to enhance our model by incorporating ecological interactions (i.e., interspecific competition), particularly for long-term predictions, as multiple dung beetle species have been introduced.Koch et al. (

1
Model parameters used to estimate dung beetle species spread and abundance.When a single value is given, it indicates that the same value was used for all species.Offspring are defined as winged dispersing adults. 1 time step (ts) = 6 months.
, FiguresS4 and S5.These parameters were used as model inputs.The spread models found the average rate of spread was ca. 30 km/year for Ont.binodis, 79 km/year for E. intermedius, 35 km/year for Oni.alexis and 28 km/year for Ont.taurus.For the average rate of spread estimation, we assumed the area occupied is a circle (see Appendix S2).The predicted numbers of individuals per grid cell 5 and 10 years after release are shown in Figure5.After 10 years, the area occupied was ca.66,874 km 2 for Ont.binodis, 99,175 km 2 for E. intermedius, 34,175 km 2 for Oni.alexis and 60,500 km 2 for Ont.taurus.See TableS8for more detail of predicted occupancy and abundance for each species.The total number of individual dung beetles predicted was ca. 12 billion Ont.binodis, 22 billion E. intermedius, 1 billion Oni.alexis and 10 billion (Ont.taurus).Given the total area of 252,098 km 2 for Sensitivity analysis of the following parameters: 'FlightDistance' (blue), 'MaxNAlates' (green), 'EggToAdultSurvival' (red) and 'ProbAdultSurvival' (orange).Graph shows ±20% from mean 'default' parameter value.anabundance of ca.47,246 beetles/km 2 for Ont.binodis, 88,100 beetles/km 2 for E. intermedius, 4,189 beetles/km 2 for Oni.alexis and 41,374 beetles/km 2 for Ont.taurus.DISCUSSIONOur model predicted the spread of four dung beetle species introduced to Western Australia.The model results highlight the variation in the predicted rate of spread among species, with the predicted average rate of spread varying from a high of ca.79 km/year for E. intermedius to a low of 28 km/year for Ont.taurus.After 10 years from the initial releases, the area occupied in the South West agricul-