Modeling round goby Neogobius melanostomus range expansion in a Canadian river system

We applied a gamma transit time model to predict the rate of range expansion of the round goby (Neogobius melanostomus Pallas, 1814) in the Trent-Severn Waterway (Ontario, Canada). Gamma distributions were fit to actual transit times of the population front from 2009 to 2011. A lack of model fit in the second year is thought to be indicative of an upstream bait bucket introduction, and this model may be useful for identifying such events. Range expansion predictions were highest in high quality habitats at 9.3 km/year, with a 5% probability that highly mobile individuals may disperse 27 km/year. The model also predicts the arrival time of the population at any distance from the population front with a given confidence interval. The estimation of a timeline for range expansion and determining underlying factors affecting the spread of invasive species could inform preventative strategies. This model is potentially useful in predicting transit times of other invasive species expanding their range in linear space, and in separating natural population expansion from additional human-assisted movement in the same system.


Introduction
Once established, invasive species often undergo range expansion as individuals disperse from the population core (Shigesada and Kawasaki 1997).This is an important aspect of species invasions because the spatial extent of their range influences the degree of their impact on recipient ecosystems (Nentwig 2007).The spread of invasive species is affected by environmental conditions and geographical limitations in relation to species biological characteristics (Shigesada and Kawasaki 1997).Understanding these characteristics aids in development of measures to prevent their spread (Sakai et al. 2001;Kolar and Lodge 2002).
Knowledge of current and potential distributions is crucial for assessing management options for invasive species (Gormley et al. 2011), but it is very difficult to determine their distribution or rate of spread because their presence often goes undetected under low-density conditions during initial stages of invasion or at the edges of their range (Vélez-Espino 2010; Jarić 2012).However, both current and future distributions can be estimated using probabilistic models (Matis et al. 1992;Sharov and Liebhold 1998;Jarić et al. 2012).The gamma transit time model was developed by Matis et al. (1992) to predict the range expansion of invasive Africanized honey bee (Apis mellifera) along the United States-Mexico border.The model uses the gamma distribution to describe observed transit times of a population front.By using transit times instead of times of transit, extreme observations of dispersal at the leading edge of population fronts can be included in the model, which can be used to make probabilistic predictions of arrival times of an invading population at locations of interest.
The round goby (Neogobius melanostomus Pallas, 1814) has been highly successful in expanding its range in both North American and European ecosystems (Jude 2001;Erős et al. 2005;Poos et al. 2009;Borcherding et al. 2011).Its spread has been attributed to a combination of human vectors (ballast water, bait buckets) and natural population expansion (Charlebois et al. 1997;Wolfe and Marsden 1998;Raby et al. 2010).Natural range expansion has been characterized as slow and continuous, while human mediated dispersal can result in colonization far from a population front (Bronnenhuber et al. 2011).However it is difficult to differentiate between natural and assisted dispersal, or to determine the natural dispersal capabilities of this species.For example, round goby range expansion was observed at a rate of 25 km/year through Chicago inland waterways from Lake Michigan to the Mississippi River (Steingraeber and Thiel 2000).In contrast, upstream range expansion averaged only 3.2 km/year in the Trent-Severn Waterway (Ontario, Canada) in the first 5 years after detection (Raby et al. 2010), and a maximum expansion rate of only 1.0 km/year was detected from Duluth Harbour, Lake Superior into the St. Louis River (Bergstrom et al. 2008).The furthest recorded dispersal distance of an individual round goby was 2 km in 218 days from a tagging study (Wolfe and Marsden 1998), suggesting that either the higher rates of observed range expansion are driven by human-mediated movement of individuals, or that longer natural migrations of individuals have not been detected by tag and recovery studies conducted to date.
In this study, we test a transit time model for predicting the rate of spread of round goby in an Ontario waterway, which incorporates the effect of habitat quality on dispersal rate.We then use the model to estimate the probability of long distance dispersal events.Finally, we examine the ability of the model to distinguish between natural and human mediated dispersal in an area upstream of an invasion front.

Study area
The Trent-Severn Waterway is a system of connected lakes, rivers and canals in a 12,550 km 2 watershed in Southeastern Ontario (Minns et al. 2004).It is a navigational waterway that contains many dams and locks, connecting Georgian Bay, Lake Huron to the Bay of Quinte, Lake Ontario.Study locations were located in Rice Lake, the Otonabee River, and Little Lake from 44.21999°N, 78.11313°W to 44.29892°N, 78.30481°W (Figure 1; for detailed site description see Brownscombe and Fox 2012).Water flows downstream from Little Lake, through the Otonabee River into Rice Lake and the Trent River.The round goby was first reported in the Trent River in 2003 at a location downstream of Lock 18 in the town of Hastings (44.31078°N, 77.952872°W) (Raby et al. 2010)

Data collection
Round goby distribution was sampled on six one-month periods between May 2009 and August 2011 (May and August of 2009 and 2010, June and August of 2011) using the same angling removal method to assess goby density as was used in 2007 and 2008.This method involves angling with small baited hooks within a 2 m 2 floating barrier for 20 minutes at each site (see Gutowsky et al. 2011).Sampling was conducted at 75 sites in 2009, and 100 sites in subsequent years at the upstream edge of its range in the Waterway, including areas beyond its detected range where future range expansion was predicted to occur.Site selection moved further upstream each year to follow the population front.Sites were randomly selected using a random point generator (geomidpoint.com).Site selection was stratified by habitat; concentrated primarily where rocky substrates were found (≈75%), since rocky substrates have been shown in previous studies to be preferred round goby habitats (Charlebois et al. 1997;Ray and Corkum 2001;Brownscombe and Fox 2012).

Range expansion model
To model the movement of the population front, a transit time (months/km) was calculated for each site where round gobies became present upstream of its detected range the year prior, calculated from the distance of the site from the previous range edge, and the date of detection.We sampled the entire expanded range for round goby presence, but our goal was to model only the pioneering population front.Therefore, we excluded sites that were clearly not part of the pioneering population front by using only sites with actual transit times of 3.5 months/km or less to predict the movement of the population front exclusively.This range of transit times was chosen based on the rates of movement of the population front observed during this study, and provided a sufficient range in values for gamma distribution fitting.
Actual transit times of round goby range expansion were categorized into bins of 1 month/km to be fitted with a gamma distribution, following the method of Matis et al. 1992.The gamma distribution can be used to fit a wide range of dispersal data (Taylor 1978;Matis et al. 1992).It is also highly tractable, and confidence intervals for arrival times can be calculated from the chi-square distribution (Matis et al. 1992).
The gamma distribution can be described by the probability density function: where β is a scale parameter with units of time, α is a unitless shape parameter, and Γ(α) is the gamma function (Johnson and Kotz 1970).This function describes the probability that x amount of time will pass before β events occur at a rate of α.As described by Matis et al. (1992), the mean E(t) and variance V(t) for transit time t were calculated as (Johnson and Kotz 1970): The model was fitted to our data (binned frequency distributions of transit times of the population front) by varying α and β parameters to achieve the same means and standard deviations (to the nearest one hundredth of a decimal place) of modeled transit times to actual transit times (Equations 2 and 3).Model probabilities were multiplied by the total frequency of actual data to obtain bin frequencies, which were compared to actual frequencies using chi square goodness of fit tests for model fit.Bin frequencies smaller than five were combined for statistical analysis, and the level of significance was p=0.05.The model was fitted to 2009−2010 data (n=22 sites at the population front) and 2010− 2011 data (n=68) independently.We also used it to determine whether a human-assisted introduction could be distinguished from unassisted range expansion on the basis of what we believe to have been a bait bucket introduction upstream of the population front.In 2009 we received a report that an angler had caught round gobies at Lock 19, 31 km upstream of the population front at the time.Long distance dispersal can occur with invasive species, resulting in the formation of separate colonies (Shigesada et al. 1995).However the natural dispersal capabilities of the round goby are relatively poor (Wolf and Marsden 1998;Hoover et al. 2003), and upstream range expansion has been consistently characterized as continuous (Bergstrom et al. 2008;Brownscombe and Fox 2012).In contrast, human assisted introductions through ballast water and bait buckets are highly prevalent in the Great Lakes Watershed (Charlebois et al. 1997;Raby et al. 2010;Bronnenhuber et al. 2011), and Lock 19 is a popular destination for recreational anglers, making a bait bucket transfer the most likely culprit.
Our sampling efforts in 2009 and 2010 did not confirm round goby presence at Lock 19, but they were detected there in our samples in 2011.Natural range expansion would have been unusually rapid from 2010 to 2011 (20.4 km compared to 9.1 km in 2009−2010), to explain the advance of the observed population front at Lock 19 in May 2011 (Figure 1).The spatial pattern of round goby density and site occupation in the area of range expansion was also atypical from previously observed population fronts.In 2009 and 2010 both site occupation and abundance decreased consistently towards the front (Brownscombe and Fox 2012), whereas in 2011, both were high at the front, followed by a large stretch of river (≈6 km) with very low site occupation and low abundance in the area of range expansion between the former front and Lock 19.Round goby were also absent from many high quality habitats in this intermediate area, which was never observed in 2009 or 2010 (Brownscombe and Fox 2012).
Round goby distribution and density were more consistent with previously observed population fronts in an area 6.5 km from the actual front (Intermediate area; Figure 1), where a gradual decrease in occupation and abundance was followed by absence from high quality habitats.As a result of these anomalies and the presumed bait bucket introduction, we generated a modified 2010−2011 distribution, in which sites above the first high quality site with gobies absent were excluded.The transit time model was then tested with 2009−2010 round goby distribution data, modified and unmodified 2010−2011 distribution data, and both years of data combined, excluding suspected bait bucket sourced sites (n=57).Finally, the dataset was separated into high quality habitat sites (>20% rock; n=40), and low quality habitat sites (≤20% rock; n=17) (Brownscombe and Fox 2012) to compare range expansion rates as a function of habitat quality.
An arrival time was predicted for a location upstream in the Waterway by multiplying α by the distance, while β remained constant (Matis et al. 1992).The 95% confidence intervals for arrival times were approximated from the 2.5 percentiles of the chi-square distribution (Matis et al. 1992).The chi-square distribution is a special case of gamma when α = ν/2 and β = 2.0, where ν represents the degrees of freedom (Johnson and Kotz 1970).The corresponding chi-square values were multiplied by β/2, to determine the upper and lower 2.5% confidence limits of arrival times.The probability of long distance dispersal events was also calculated from the cumulative probability of the model.

Data analysis
Round goby transit times included in the model were compared between time periods (2009− 2010, 2010−2011), inclusion and exclusion of bait bucket introduction sites (2010−2011 only), and high quality vs poor quality habitats assessed over the 2009−2011 period.Comparisons were made with one-way Analysis of Variance and Tukey's HSD post-hoc analysis.Data were log10-transformed prior to analysis, and the level of significance was set at p < 0.05.

Results
The model was a good fit to the 2009−2010 range expansion data, but was a poor fit to the 2010−2011 unmodified dataset (Figure 2a,b;  Table 1).In contrast, by excluding round goby occupation sites suspected to have occurred from an upstream bait bucket introduction, the model fit the 2010−2011 data well (Figure 2c; Table 1).
It also fit the two-year combined dataset when excluding these sites (Figure 2d), and when separating this dataset into high quality and low quality habitat types (Figure 2e,f; Table 1).
Actual range expansion was more rapid from 2010 to 2011 (12.9 km/year excluding hypothesized bait bucket sourced sites) than from 2009 to 2010 (10.4 km/year) (Figure 1), which resulted in a strong positively skewed transit time probability distribution in the second year (Figure 2).Transit times varied significantly between years (F 5,233 = 18.2, P <0.001) and were significantly faster in 2010−2011 than in 2009−2010 (Tukey's HSD, p < 0.05).The 2010−2011 model predicts a 52% faster rate of spread than the 2009−2010 (Table 1).By combining these datasets and modeling range expansion from 2009 to 2011 (Figure 2d), the arrival of round gobies (in detectable numbers) at Lock 22 (8.1 km upstream in the Waterway from the round goby population front in August 2011) is predicted to occur in 1 year, but may occur as soon as 8.4 months (0.7 years) based on the 95% confidence interval (Table 1).The use of only higher quality sites in the model predicts a 24% faster rate of spread than that predicted when only the lower quality habitat sites were used, and the former predicts 9.3 km/year of range expansion (Table 1).The high quality habitat model predicts arrival at Lock 22 as soon as 7.2 months (0.6 years) (Table 1).By examining the cumulative probability distribution for dispersal into high quality habitats (the most rapid range expansion), the model predicts that there is a 95% probability that at least some individuals will disperse 4.9 km/year upstream in the Waterway, while there is a 5% probability that some highly dispersive individuals may travel 27 km/year (Figure 3).

Discussion
The gamma transit time model was shown to be a useful tool for predicting many aspects of round goby range expansion in the Trent-Severn Waterway.An important benefit of this model is that it predicts the arrival of round gobies at specific locations based on their distance from the population front; a potentially useful tool for informing where management should be directed.There are many types of physical and nonphysical barriers (e.g., electric barriers, velocity barriers, acoustic deterrents) that have been used in the past to reduce the spread of a variety of pest species (Dahlsten and Garcia 1989;Liebhold et al. 1992;Noatch and Suski 2012).For obvious reasons, barriers must be located beyond the expanding population front (Sharov and Liebhold 1998), so it is critical to develop methods of predicting its location.These fronts are often populated by a few dispersers below the minimum detectable density, but models can be used to identify their probabilistic front locations.
Habitat can be an important variable affecting the spread of invasive species (Shigesada and Kawasaki 1997).Round gobies have an affinity for rocky substrates (Brownscombe and Fox 2012), and the model applied here predicts more rapid range expansion into these high quality habitats.Therefore population spread would be expected to occur faster in ecosystems with a high proportion of rocky habitats.Large areas of poor quality habitat have been found resistant to round goby invasion, with much slower rates of colonization (Cooper et al. 2007;Young et al. 2010), which is consistent with our model predictions.
The gamma transit time model also provides the utility of predicting long distance dispersal events, which is particularly useful for invasive species such as the round goby where the probability of establishment is high even under lowdensity conditions (Vélez-Espino et al. 2010).However, because the probability distribution used to predict dispersal distances is asymptotic, it can produce unrealistic estimates, and the energetic limitations to movement of the species must be considered.Large round gobies have been found to exhibit swimming speeds of up to 20 cm/s at 17°C., but only for a relatively short duration of time (up to 72 minutes; Hoover et al. 2003).Hypothetically, if a round goby could disperse at this speed for 12 hours per day, this individual could travel up to 26 km in a month.However, little is known about the dispersal capabilities of highly mobile individuals, and this is an important parameter for predicting population expansion.Regardless, the tail of the probability distribution provides some insight into the scale and probability of dispersal distances for this species, and management strategies should begin much sooner than the predicted arrival of a detectable population front.
Sampling low-density populations is difficult due to the low probability of catching the targeted species.By using the gamma distribution to model transit times, we were able to identify patterns consistent with low-density bait bucket introductions that would normally escape notice.In our study, the poor fit of the gamma distribution model to unmodified 2010 to 2011 range expansion data was caused by an atypical distribution at the population front.Because site occupancy was higher at the pioneering population front than in the intermediate area of range expansion, the skew of the transit time distribution was too strong to allow a good fit with the model (Figure 2b, Table 1).Once those sites estimated to be sourced from an upstream introduction were removed, the model fit the distribution well (Figure 2c).These sites were removed from the dataset using a criteria based upon the known spatial ecology of this species at population fronts (see Methods for description), not to fit the model.Therefore our results suggest that this model may be useful for detecting human-assisted introductions when the spatial separation between the introduction site and infested zone makes the introduction uncertain.In this study, we were able to identify a probable human-assisted dispersal event through extensive sampling.However, this is not always the case with invasive species, and this model may be useful in these instances.By distinguishing between natural and assisted dispersal, we can gain a better understanding of how invasive species are spreading and improve management strategies.
Invasive species are known to expand their ranges by both short distance dispersal into neighbouring habitats, and long distance dispersal, which can result in the formation of colonies that are spatially separated from the previous population front (Shigesada et al. 1995).The latter may result in spatial distributions similar to human assisted introductions, which may also cause a lack of fit to a gamma transit time model.In the case of the round goby, colony formation is more likely to occur in downstream range expansion due to larval drift, whereas upstream range expansion has been consistently characterized as continuous (Hensler and Jude 2007;Bergstrom et al. 2008;Brownscombe and Fox 2012).Based on knowledge of upstream range expansion mechanisms and dispersal capabilities of the round goby (Wolf and Marsden 1998;Hoover et al. 2003), the formation of a colony 31 km upstream of the population front in the Trent-Severn Waterway is far less likely than a human assisted introduction.However, the potential for colony formation should be strongly considered when attempting to separate natural from assisted dispersal with a gamma transit time model.Assisted dispersal has been suspected as a strong contributor to secondary round goby range expansion throughout the Great Lakes Basin through circumstantial and genetic evidence (Corkum et al. 2004;Bronnenhuber et al. 2011).This combination of natural and human-assisted dispersal is defined as stratified dispersal, which is recognized as an important component to range expansion of many invasive species (Parisod and Bonvin 2008;Darling and Folino-Rorem 2009;Bjorklund and Almqvist 2010).
The rates of range expansion of invasive species are often highly variable, and can be influenced by environmental characteristics, propagule pressure, and population dynamics (Shigesada and Kawasaki 1997;Lockwood et al. 2005).The density of an invasive species may be an important driver of range expansion, particularly for the round goby, which has very high site fidelity (Wolf and Marsden 1998).As the population size becomes larger, and local habitats become saturated, range expansion may become more rapid.Indeed, round goby range expansion started slowly and accelerated over time in the Trent-Severn Waterway since its introduction in 2003, but has become more consistent over the last 3 years (Brownscombe and Fox 2012).Therefore it is important to consider the population density and the invasion timeline when making predictions with a gamma transit time model, as it does not include density as a factor, but makes predictions based on observations of advancements of the population front.However, the simplicity of the model applied here has certain advantages for future applications.Reliable density data is often difficult to acquire, particularly for the round goby, where estimates vary greatly from different sampling methods (Johnson et al. 2005).This model only requires observations of advances of a population front.
In conclusion, the gamma transit time model is a useful tool for predicting many aspects of round goby range expansion in a linear system, including the prediction of arrival times at upstream locations of interest, and the probabilistic prediction of long distance dispersal events.Our test of the model suggests that it is also useful in distinguishing between natural and humanassisted range expansion, but further testing is required using experimental introductions to determine the conditions under which humanassisted introductions can be detected.

Figure 1 .
Figure 1.Round goby distribution in a section of the Trent Severn Waterway from 2003 to 2011.Black hatched lines indicate extent of upstream distribution in a given year, black circles are navigational locks, and stars are points of introduction.Sampling areas in 2009/2010 are shown in light grey, 2011 in dark grey.

Figure 2 :
Figure 2: Probability distribution of actual (bars) and model (line) transit times of round goby range expansion from (A) 2009 to 2010, (B) 2010 to 2011 including sites suspected to be sourced from an outside introduction 31 km upstream of the infested zone in 2009, (C) 2010 to 2011 excluding outside sourced sites, (D) 2009 to 2011 excluding outside sourced sites, (E) 2009 to 2011 low quality habitat sites, and (F) 2009 to 2011 high quality sites in the Trent-Severn Waterway.Hatched line indicates poor model fit.

Figure 3 .
Figure 3. Probability of round goby dispersal distances (km) from the population front predicted by the 2009 to 2011 high quality habitat model of round goby range expansion in the Trent-Severn Waterway.Dashed line indicates 5% probability.

Table 1 .
Actual transit times for round goby range expansion, model parameters, fit, predicted arrival times (years) of a detectable round goby population at an upstream location at Lock 22 (95% CI), and rate of spread (95% CI).Italics indicate predictions from a model with poor fit to actual data.