Predicting the outcomes of management strategies for controlling invasive river fishes using individual- based models

1. The effects of biological invasions on native biodiversity have resulted in a range of policy and management initiatives to minimize their impacts. Although management options for invasive species include eradication and population control, empirical knowledge is limited on how different management strategies affect invasion outcomes. 2. An individual- based model (IBM) was developed to predict how different removal (‘culling’) strategies affected the abundance and spatial distribution of a virtual, small- bodied, r - selected alien fish (based on bitterling, Rhodeus sericeus ) across three types of virtual river catchments (low/intermediate/high branching tributary configurations). It was then applied to nine virtual species of varying life- history traits ( r- to K- selected) and dispersal abilities (slow/intermediate/fast)


| INTRODUC TI ON
Biological invasions are an important component of global change (Simberloff et al., 2013), which have resulted in the application of a range of policy and management initiatives to minimize their impacts (Larson et al., 2011). Optimising population control efforts within management programs remains highly challenging, and considerable uncertainty remains in how to apply limited staff and financial resources to population control measures, especially as to how these measures should be applied spatially and temporally (Maguire, 2004).
Following the introduction of a new species, the probability of an invasive population developing depends on the interactions of a range of biotic and abiotic factors, including the species' dispersal abilities and life-history traits, and the environmental conditions encountered (Catford et al., 2009). Dispersal rates are important and are influenced by habitat connectivity (Hastings et al., 2005), which is often more constrained in freshwaters than in terrestrial environments (Gozlan et al., 2010;Radinger et al., 2017). Management responses also influence invasion probabilities , and their effectiveness generally increases when they occur soon after introduction, when the species is spatially constrained and of relatively low abundance (Rytwinski et al., 2019).
Effective strategies for controlling invasive species must consider two important aspects. First, the most efficient, cost-effective and safe means of removing individuals from a local population must be determined, which will vary by the invader's life stage and the ambient environmental conditions (Buhle et al., 2005). Second, the timing and location of the control measure needs determining. For example, in invasive plant management, there is debate about whether efforts should be focused at the range front, where the invader is of low abundance and patchy in distribution, or at the invasion core, where populations are established and usually more abundant (e.g. Hastings et al., 2006;Januchowski-Hartley et al., 2011). Invader dispersal rates and abundances can also be affected by habitat complexity, with terrestrial invaders generally spreading more easily in unfragmented than fragmented landscapes (Dewhirst & Lutscher, 2009). Therefore, understanding the ease with which invaders spread in different environmental configurations needs consideration within management planning (Lurgi et al., 2016).
In managing freshwater invaders, the utility of individual-based models (IBMs) to predict rates of establishment and spread has recently been highlighted (Dominguez Almela et al., 2020), along with comparisons of the methods for controlling and/or eradicating local populations (Rytwinski et al., 2019). There is, however, less known on how invader dispersal rates affect the long-term efficacy of population control, especially at large spatial scales and in open, linear systems. In rivers, while controlling invasive fish can be a management priority, their removal usually relies on capture methods (e.g. nets, traps, electro-fishing; Rytwinski et al., 2019). Although these are effective at capturing (and culling) fish species across most of their size range (Davies & Britton, 2015), they are unlikely to remove all the population, and surviving individuals can then potentially compensate for losses (Berry et al., 2012). This can result in culled populations rapidly recovering to their previous abundances and management objectives not being met (Davies & Britton, 2015;Dominguez Almela et al., 2020).
Here, our aim was to develop an IBM to predict the effects of management control strategies on invasive river fishes. We applied it within factorial experiments to 'virtual' invasive fish with a range of life-history traits in river networks of varying spatial complexity to predict how different culling strategies affected fish spread and abundance. The objectives were to predict: (a) how altering the application of culling (rate/location/timing/life stage) affected the abundance and spatial distribution of an initial invasive fish, based on the demographic characteristics of bitterling Rhodeus sericeus, a small-bodied invasive fish (Dominguez Almela et al., 2020); (b) how invader abundances and distributions varied across nine virtual invasive fish differing in life-history traits and dispersal abilities; and (c) the trade-offs between management effort (cull rate and the number of patches culled) and the predicted invader abundances and distributions to optimize management responses. The final model was then applied to a real-world scenario to predict how management interventions could have constrained the real-world invasion of bitterling in the River Great Ouse, Eastern England (Dominguez Almela et al., 2020).

| Model configuration and virtual species
The model was implemented in a customized version of the individualbased spatially explicit modelling platform RangeShifter (Bocedi et al., 2014(Bocedi et al., , 2020, incorporating a new module for managing invasive 5. Synthesis and application. This IBM predicted how management efforts can be optimized against invasive fishes, providing a strong complement to risk assessments.
We demonstrated that for a range of species' characteristics, culling can control and even eradicate invasive fish, but only if consistent and relatively high effort is applied.

K E Y W O R D S
biological invasion, dispersal, RangeShifter, river catchment, simulation model species. River catchments of an overall mean extent of 74.5 ha and of three basic types were created using ArcGIS Pro, where the configuration varied by the extent of its tributary branching: high branching (12 tributaries), moderate branching (eight tributaries) and low branching (four tributaries). Each catchment type was replicated three times in slightly different configurations ( Figure 1). All catchments were created in raster format at 50-m resolution, and each cell was assigned a habitat quality score proportional to stream width (Dominguez Almela et al., 2020). Catchments comprised a mean of 298 (SD: ±20) cells divided into a similar number of non-overlapping patches (mean = 44.9 patches/catchment, SD = ±0.3), which were delineated manually such that no patch overlapped a confluence. The ratio of good to poor quality habitat sections was similar in all configurations; in the main stem, it was kept consistent at 1:1, and in the tributaries, the ratio was maintained across replicates within the branching groups (low/moderate/ high) at an overall mean of 8:15.
The initial virtual fish, whose demographic and dispersal parameters were drawn from bitterling (cf. Dominguez Almela et al., 2020), had the characteristics of an alien fish that followed a stage-structured population dynamic: juveniles (<1 year old), subadults (1-2 years) and adults (over 2 years). Its population parameters (Table 2; Table S1) were also similar to other small-bodied, invasive, pest fish species, such as topmouth gudgeon Pseudorasbora parva (Britton & Gozlan, 2013). It was set to reproduce once per year, exhibit density-dependent fecundity and have a sex ratio at birth of 1M:1F. Its survival probability was also density dependent, applied and weighted per stage, so the effect of subadults on the survival of adults was 10% of the effect of adults on each other and on subadults (Dominguez Almela et al., 2020). After reproduction, juveniles, subadults and adults could disperse according to a density-dependent emigration probability, limited to one dispersal event per lifetime. Movement from the natal patch during the transfer phase of dispersal was modelled by the stochastic movement simulator (SMS; Palmer et al., 2011), which simulates movement from cell to cell on the basis of perceived costs within a limited perceptual range and a tendency to follow a correlated path (directional persistence).
A relative cost of movement map was derived from the habitat map, such that perceived costs were inversely related to habitat quality; thus, upon reaching a confluence, a disperser would more likely move into the wider of the two streams available. Individuals could settle in any non-natal patch subject to an inverse density-dependent settlement probability ( Table 2; Table S1).
For all model simulations, initial populations were established in 10 patches (including in tributaries) at the upstream end of the catchment (Figure 1). These were mainly in first and second order streams, the rationale being that this enabled the invasion front to spread mainly downstream into streams of higher order (Kim et al., 2021), although the direction of individual movement was stochastic in the model and could occur in either direction unless the fish was in a terminal patch. Propagule pressure by the random addition of further individuals into the system was not considered in the model. Each simulation (a single combination of culling parameters and catchment; Table 1) was run for 30 years and replicated five times. A 'control' simulation was run for each catchment in which no management was applied.

| Experimental designs
The management of the modelled populations took the form of an annual cull, the format of which was controlled by a series of variable F I G U R E 1 Virtual catchments having a low number of tributaries (a1-3), medium (b1-3) and high (c1-3). Grey thick lines highlight the initialized patches in the upper areas of the catchment, and the arrows indicate the direction of flow along the main stem parameters. A specified number of patches was selected each year for culling in which individuals were subjected independently to a mortality probability determined as a stage-dependent logistic function of density and subject to a maximum culling rate. This means that the proportion of individuals removed from a patch per year (the 'cull rate') would depend on its density, and low invader numbers would translate into low removals, reflecting the difficulties of catching those individuals when the population is very low. In contrast, if the density was high, the logistic function would result in removal of a higher proportion of individuals up to a pre-defined maximum. Both the culling rate and the number of patches selected were kept constant across years. We modelled the cull as occurring after dispersal and its spatial application as biased towards recently colonized patches based on results from preliminary experiments (see Appendix S1-S3). Then, the first experiment assessed the interaction of culling rate and number of patches culled by varying the cull rate on an initial virtual species across different numbers of patches (Experiment 1, Table 1). The second experiment varied the demographic (as type of demographic species, SpType) and dispersal (as type of dispersal species, SpDispType) characteristics of the model alien species that was being managed (Experiment 2, Table 1).
Finally, in the third experiment, we applied the model to a real case study using bitterling as the model species and the Great Ouse as the model catchment, following the study by Dominguez Almela et al. (2020).

| Experiment 1: Interaction of culling and management effort
Six levels of the maximum number of patches culled (N), ranging from low to high numbers (4 to 24) of patches culled per year, were tested against eight levels of maximum cull rate (CR, 0.2-0.9; Table 1). The cull was applied to all three life stages simultaneously. The experiment was run for all nine catchments (432 parameter combinations).
Following the model simulations, the results of Experiment 1 were analysed using two population-level summary statistics that were extracted from the model output data: where Nind t is the total number of individuals across the catchment at time t and n is the number of years during which there was active growth. The simulated population when no culling was applied was used to determine n, the number of years of approximately constant growth before the population growth trajectory began to decline.
where NOccupPatches t is the number of occupied patches at time t and n is the number of years for the period of active growth as above.
2.2.2 | Experiment 2: Trade-off between the number of patches culled and the cull rate within each patch across a range of species This experiment evaluated a potential resource-limited trade-off between the maximum number of patches culled and the maximum cull rate within each patch across a further eight contrasting virtual species (specific culling strategy; SCS), to identify whether patterns detected in the initial species were common across species with contrasting life-history traits and dispersal abilities (Table 1). It was based on the assumption that the total management resource (finance, manpower) was fixed annually, and so managing more patches in a year would result in reduced effort per patch, reducing the effects of the culling on the invader. For example, culling all patches would spread resources thinly across the whole catchment, so the lowest maximum cull rate (CR = 0.1) would have to be applied per patch in each year in this scenario (Table 1). The initial virtual species (Experiment 1) was the reference species, having intermediate demographic traits and dispersal abilities, and eight additional virtual species were developed that varied by their: (a) demographic traits and/or (b) dispersal traits.
(a) ranged from being more intensely r-selected (i.e. higher fecundity, lower survival) to being more intensely K-selected (lower fecundity, higher survival). (b) Varied from being better dispersers by increasing stage-dependent maximum emigration probabilities (more dispersers per generation) and decreasing maximum settlement probability (more likely to keep moving), to poorer dispersers by decreasing stage-dependent maximum emigration probabilities (fewer dispersers per generation) and increasing maximum settlement probability (less likely to keep moving; Table 2). The experiment was applied to all nine catchments (729 parameter combinations).
As there was substantial variation between population trajectories in the experimental predictions, fixing a specific year as the basis for calculating the rate of population increase (P 1 ) and change in patch occupancy (Q 1 ) was inappropriate. Therefore, for each replicate simulation trajectory, the rate of population increase in the first decade (P10) was calculated, using: where pop[year=10] is the number of individuals at year 10 and pop [-year=1] is the number of individuals at year 1, and similarly for years 6-15 inclusive to give P15; 11-20 (P20), 16-25 (P25) and 21-30 (P30).
Finally, the maximum mean annual increase in population size achieved was determined as: The same approach was used to calculate the maximum decadal rate of change in patch occupancy (Q 2 ): (1) Rate of population increase P 1 :

| Statistical analyses
For both P n and Q n , factorial linear models for each experiment were fitted using R version 3.6.3 (R Core Team, 2020) to partition the variance in the response variable. The models incorporated, as appropriate, the management scenarios that were applied (N, CR and SCS), the three sets of river catchments (low, medium and high branching; B, ID) and species (SpID, SpType and SpDispType) ( Table 1).
The effect of the factor(s) that had the greatest influence on model outcomes was investigated further using posterior marginal means analyses (package 'emmeans', Lenth, 2020).

| Experiment 1: Interaction of culling and management effort
The predicted effects of culling on both P 1 and Q 1 increased as both the maximum cull rate (CR) and maximum number of patches (N) increased, which together accounted for most of the variation (Table 3; Tables S2b and S3b). When CR was low, it had only minor effects on P 1 and Q 1 , except when N was high ( Figure 2). When CR was increased to a medium level (0.5, 0.6), then values of P 1 and Q 1 suggested that culling could contain the spread of the invader, but when CR was ≥0.7, eradication of the species was possible if N was also high (at least 16 patches culled per year; Figure 2). All replicates led to eradication when CR was at least 0.8 and N was ≥12 patches per year, and also when CR was 0.9 and N = 8 patches per year.

| Experiment 2: Trade-off between the number of patches culled and the cull rate within each patch across a range of species
In general, the major source of variation in P 2 and Q 2 was from differences between the characteristics of the nine species (SpType, SpDispType and SpID; Table 4; Tables S2b and S3b). Marginal means analysis between species type (SpType and SpID) and the SCS revealed that P 2 and Q 2 could be substantially reduced when N was set at a relatively high number of patches, despite culling being at less than maximum efficiency (Figure 3; Figure S4). Five out of the nine species presented at least two scenarios in which the virtual species was eradicated in some of the replicates (  (Table 5). In particular, the scenario where 0.8 cull rate was applied to nine patches predicted that for the r-selected slow dispersers, there were no replicates where the species was not eradicated. In r-selected fast dispersers, P 2 and Q 2 were reduced, but their populations were only predicted to be eradicated in a single replicate (Figure 3; Figure S4). Conversely, 18% of intermediate and 61% of slow dispersers within K-selected species were successfully eradicated in simulations at CR ≥ 0.7, while fast disperser abundances were reduced, but were never fully eradicated at any of the nine CR scenarios (Table 5; Figure 3; Figure S4).
There were also marked differences in the response of P 2 and Q 2 in relation to the different dispersal abilities of the nine species (SpDispType; Figure 4). Species with slow dispersal abilities required less effort to control, enabling lower numbers of patches to be culled to achieve similar outcomes as intermediate or fast dispersers TA B L E 2 The model parameters of the nine virtual species used in Experiment 2. See Table S1 for parameters common to all species Max. emigration probability in adults 0.14 0.14 0.14 Max. settlement probability 1 1 1 Per-step mortality 0.01 0.01 0.01 *Species 5 was the same as the single species used in Experiment 1 (and in preliminary experiments 0a, 0b in Supporting Information).; **1/b is the rate of density dependence, that is the rate at which mean fecundity decreases with increasing local density.
at higher numbers of patches. Moreover, the intermediate and fast dispersers revealed some compensatory responses to the culling, suggesting their invasion could actually benefit from some culling scenarios through increased P 2 and Q 2 .

| Experiment 3: Case study on specific management to control bitterling in the River Great Ouse
Our IBM predicted a significant response in this bitterling population to the different culling scenarios in both P 2 (F 8,2,241 = 165.23, p < 0.001) and Q 2 (F 8,2,241 = 1,953.5, p < 0.001; Figure 5) analogous to that predicted in Experiment 2 for intermediate demography/dispersers (cf. Figure 4). Where Dominguez Almela et al. (2020) predicted bitterling would occupy 90% of patches in the river in 2045 (Figure 6a), the application of a yearly SCS of CR = 0.7 and N = 14 (5% of patches) was predicted here to result in the population occupying only 21% of the area, which was similar to their spatial extent recorded in 1984 (Figure 6b).

TA B L E 3
Principal sources of variance explained (%) in the summary statistics rate of population increase (P 1 ) and rate of change in patch occupancy (Q 1 ) for Experiment 1. Factors: maximum number of patches culled (N), maximum cull rate (CR), catchment branching (B) and catchment ID number (ID)

F I G U R E 2
Interaction effects of the maximum cull rate (CR) and maximum number of patches culled (N) on the rate of population growth (P 1 ) and change in patch occupancy (Q 1 ) in Experiment 1

| D ISCUSS I ON
Our model predictions provide a series of novel insights that should assist decision-making for managers dealing with invasive fish specifically and invasive species more generally. Our IBM identified the optimal management effort (as specific culling strategies involving the cull rate and its spatial application) for invaders across a range The predictions indicated that a reduced cull rate, even when applied over larger spatial areas, was generally ineffective in constraining the fish dispersal and population growth rates. Correspondingly, the application of low cull rates across large spatial areas is not an effective management option. In contrast, relatively high cull rates often predicted population eradication, even when the number of patches culled were not necessarily high (e.g. a cull rate of 0.9 in only eight patches). These predictions are thus important in the context of management planning of how specific culling strategies could be applied spatially and in relation to the ability of capture methods to remove high proportions of the target species within managed patches. They also demonstrated the compensatory responses that will occur within the fish population when the cull rate is too low, whereby the reduction in population size results in subsequently higher reproductive rates and abundances, and then faster dispersal (Berry et al., 2012).
Although the IBM was effective at predicting how the cull rate affects the abundance and dispersal of the invader, determining the extent of a fish population that can be removed effectively from a population using capture methods, such as electric fishing, can be difficult. While electric fishing is considered an effective fish capture technique, it has inherent issues relating to species detectability and detection bias (Beaumont, 2016). Its probability of capture is species dependent with, for example, it varying between 0.35 and 0.64 across 15 stream fish species (Reid et al., 2009). It can also be relatively ineffective at capturing early life stages of fish, where alternative methods might be more effective, such as micromesh seine nets and/or traps (Britton, Pegg, et al., 2011;Nunn et al., 2001).
Nevertheless, electric fishing is regularly and successfully used to remove alien fish from invaded waters when chemical treatments cannot be applied, with the meta-analysis of Rytwinski et al. (2019) reporting a 58% success rate for population eradication and 56% for population control over a range of alien fish species.
The predicted outcomes of the specific culling strategies on the invading populations were also strongly influenced by the invader's demographic traits. The initial virtual species, based on bitterling, had a suite of demographic traits that were similar to those of small-bodied invasive fish more generally that are considered pests in many parts of the world and so often receive considerable management attention (e.g. Britton & Brazier, 2006;Britton et al., 2010).
However, invasive fish with demographic traits that are less intensively r-selected or even K-selected are also considered undesirable in many parts for the world. For example, relatively large-bodied alien fish of the Salmonidae family are often targeted for management in North America (Rytwinski et al., 2019) and common carp Cyprinus carpio have also received considerable management attention in many areas of the world, including Australia (Pinto et al., 2005) and South Africa (Davies et al., 2020). It was thus important to understand how these initial predictions varied according to demographic traits and dispersal abilities. The predictions revealed that population eradication was achievable in five of the nine virtual species, including species with both r-and K-selected traits, but only when the cull rate was relatively high (≥0.7) and when at least 14 patches were culled. However, cull rates below 0.7 never resulted in eradication.
That invasive fish with r-selected demographic traits were predicted to be the most challenging to control and eradicate is arguably an intuitive result, as these traits, common in many invasive fish, enable the rapid development of highly abundant populations where individuals can then disperse (Gozlan et al., 2010). The eradication of two strongly K-selected virtual species was predicted, but only when they were intermediate or slow dispersers and at a high cull rate. The real-world simulation of bitterling in the River Great Ouse indicated that although population eradication would not occur under any of the specific culling strategies that were simulated, strong containment was possible.
These invasion management predictions should be applied to informing real-world invasion risk assessment processes (Vilizzi et al., 2019), which are fundamentally important in prioritizing the species and habitats for management (Roy et al., 2018). Best practice guidance on risk assessments for alien species indicates that two of the four main assessment components are their probable spread and impact (Roy et al., 2018). In our summary analysis of the IBM output data, we directly predicted invader spread, and the parameter P n can potentially predict impact due to the relationships of invader abundance with impact (Jackson et al., 2015). Consequently, provided some basic knowledge of the demographic and dispersal traits is available for a specific species, our IBM predictions should inform the risk assessment responses on their spread and impact, especially where there is strong understanding of the invader's demographic traits and dispersal abilities. This is because there were some strong differences in the responses of r-and K-selected species, and across varying dispersal abilities, to the different culling strategies. For example, 82% of slow-dispersing r-selected species were eradicated at cull rates above 0.7% versus 61% of K-selected species, while eradication of fast-dispersing species was generally predicted to be unlikely. Although detailed information on the dispersal abilities of alien fish is often lacking, dispersal patterns are also likely to be catchment specific, varying according to, for example, the flow regime, river network complexity and habitat suitability, and the extent of human modification of them (Caiola et al., 2014).
Correspondingly, following the detection of a new alien fish within a river catchment, practitioners should rapidly assess its life-history traits (either directly or through literature review), its ability to disperse in the catchment (in relation to both the species' dispersal characteristics and the complexity of the river network) and its current spatial extent . Concomitantly, the resource available for the control effort needs to be quantified, with identification of the desired management outcome (eradication vs. control).
In combination with the predictions outlined here across the different fish life histories and dispersal abilities, and the culling gradients, this information should then enable more informed decision-making on the actual strategy to be implemented.  (Kim et al., 2021). In reality, the release points of alien species are more stochastic, such that some dispersing aliens move from downstream to upstream areas (Vitule et al., 2012). The model also prohibited an individual from dispersing more than once in its life (Bocedi et al., 2014), an assumption potentially violated by some riverine fishes (Fausch et al., 2002;Radinger & Wolter, 2014). Moreover, the model was based on a discrete introduction of fish, and while invasions can occur from a single release event, multiple releases can also occur, which can increase the probability of invasion success as it overcomes issues such as founder effects (Lockwood et al., 2005). However, a key component of eradication attempts of invasive species is preventing their reintroduction into the treated area, and so where management efforts are ongoing to control invaders, these efforts should also include increased regulation and surveillance that aim to prevent further releases . Finally, throughout the model, the total specific culling strategy was fixed across all years, as this enabled the model to account consistently for the relationship between cull rate and the number of patches culled. However, it is acknowledged that in reality, resource availability might vary by year and, for high priority species, be increased in the short term. Moreover, it has already been discussed that the actual effort required to achieve F I G U R E 5 Effects of specified culling strategy on the rate of population increase (P 2 ) and change in patches occupancy (Q 2 ) of the bitterling species. Grey bars are 95% confidence intervals for the marginal means the higher cull rates might exceed the levels of effort used in the model due to the increased degree of difficulty of capturing fish from small populations (Britton, Pegg, et al., 2011). Consequently, while there is high confidence in our predictions of how the different cull rates applied across varying number of patches culled will affect the population of the target species, it is suggested that some caution is applied when considering that these strategies will all be of equal cost. Despite these issues, we argue strongly that the strength of our model is its ability to simulate the outcome of management interventions to control populations of invasive species that cannot be provided empirically (Dominguez Almela et al., 2020). It thus represents a major step forward in understanding how to develop strategic approaches for managing alien species in the environment.
In summary, this work provided insights into the outcomes of different control efforts on invasive fishes, and highlighted that, depending on the species characteristics, and the specific culling strategy, these outcomes can vary in target populations, but can include eradication when the target species is of low or intermediate dispersal ability and when a high cull rate is applied. These predictions strongly complement existing invasion risk assessments, and demonstrate that individual-based models are powerful tools for predicting optimal management interventions for high-risk invaders.

ACK N OWLED G EM ENTS
The RangeShifter software and manual can be downloaded from

CO N FLI C T O F I NTE R E S T
None of the authors have any conflict of interest. J. Robert Britton is an Associate Editor of Journal of Applied Ecology but took no part in the peer review and decision-making processes for this paper.

AUTH O R S ' CO NTR I B UTI O N S
All authors contributed to the study design. V.D.A. led analyses and writing, assisted by S.C.F.P., P.K.G., D.A., J.M.J.T. and J.R.B. All