Invasion of American bullfrogs along the Yellowstone River

The American bullfrog (Lithobates catesbeianus) is a globally distributed invasive species that was introduced to the Yellowstone River floodplain of Montana. Knowledge about floodplain habitat features that allow for bullfrog persistence and spread will help identify effective control strategies. We used field surveys in 2010, 2012 and 2013 to describe bullfrog spread in the Yellowstone River floodplain and the habitat features that are associated with bullfrog occupancy and colonization. Bullfrogs in our study area expanded from ~ 60 km in 2010 to 106 km in 2013, and are spreading to upand downstream habitats. The number of breeding sites (i.e., presence of bullfrog eggs or larvae) increased from 12 sites in 2010 to 45 sites in 2013. We found that bullfrogs were associated with deeper waters, emergent vegetation and public-access sites, which are habitat features that characterize permanent waters and describe human-mediated introductions. Control strategies that reduce the hydroperiod of breeding sites may help to limit bullfrog persistence and spread, while an increase in public outreach and education may help prevent further bullfrog introductions at public-access sites.


Introduction
Introduced American bullfrogs (Lithobates catesbeianus [Shaw, 1802]; hereafter, bullfrog) have been implicated in the declines of multiple amphibian and reptile species around the globe (Ficetola et al. 2007).Their large size, high mobility, generalist eating habits, high fecundity and function as a disease vector makes the bullfrog an extremely successful invader and a threat to biodiversity (Lowe et al. 2000;Nentwig 2007).An improved understanding of bullfrog invasion dynamics is needed to limit their spread and minimize impacts on native taxa (Ficetola et al. 2007).
Bullfrogs are native to eastern North America (Bury and Whelan 1984), but they now occur in most western states and provinces in North America.The impetus for bullfrog introduction seems largely to be for recreational hunting and human consumption (Bury and Whelan 1984;Jennings and Hayes 1985).Once introduced, suitable habitat in the form of permanent waters is a critical factor in bullfrog establishment and spread (Ficetola et al. 2007;Fuller et al. 2011).Bullfrogs in much of western North America require permanent water for reproduction because larvae generally overwinter at least once before metamorphosis (Cook et al. 2013).
Bullfrogs are a relatively new invader in the Yellowstone River floodplain of Montana (MT), USA.They were first documented in a permanent pond adjacent to a Yellowstone River tributary near the city of Billings in 1999 (Figure 1; Observation IDs 10028191 and 10028199 in Montana Natural Heritage Program 2014).In 2005, bullfrogs were observed again at this permanent pond but also 10-km downstream in a private campground pond and 40-km downstream in an off-channel canal next to a popular public angling site (Observation IDs 10306890 and 10305877 in Montana Natural Heritage Program 2014; Figure 2).The potential for bullfrog spread is high because this floodplain is a mosaic of human-modified permanent waters, canals, backwaters and side-channels (e.g., Fuller et al. 2011).Post-metamorphic bullfrogs can disperse long distances and are adept at colonizing new water bodies (> 1200 m; Willis et al. 1956).The Yellowstone River floodplain provides habitat for native amphibians that can be impacted by bullfrogs; including the Northern leopard frog (Lithobates pipiens [Schreber, 1782]), Woodhouse's toad (Anaxyrus woodhousii [Girard, 1854]) and Great Plains toad (A. cognatus [James, 1823]).Bullfrogs have been implicated in amphibian declines in other areas, including Northern leopard frogs (Hammerson 1982;Hayes and Jennings 1986;Johnson et al. 2011).
Here, we describe bullfrog expansion and occupancy in the Yellowstone River floodplain.We assessed bullfrog expansion with field surveys in 2010 and resampling in 2012 and 2013.To describe occupancy, we related bullfrog occurrence from the 2012 and 2013 surveys to habitat features related to human-mediated introductions, humanmodified habitats and habitat quality.

Study area
The Yellowstone River emerges from Yellowstone National Park in southwestern MT, flows through MT's northern Great Plains and drains into the Missouri River in western North Dakota (Figure 1).Much of the aquatic habitat in and adjacent to the floodplain has been modified by irrigation withdrawals, agriculture, channelization and development.In 2010, we conducted baseline surveys along the Yellowstone River from downstream of the town of Laurel (river km 611) to Finch (km 394).In 2012 and 2013, we surveyed from downstream of Park City (km 619) to Custer (river km 482; Figure 1).Surveys were conducted in the Yellowstone River, side channels, backwaters and off-channel impoundments.

Data collection
2010.We conducted visual encounter surveys (VES) and calling surveys between July 12 th and October 6 th at 110 sites.We visited ten sites 2 times per month for 3 months and all other sites once.These sites were selected non-randomly and largely based on public accessibility and prior bullfrog sightings.We focused on habitats that had structure (aquatic vegetation or woody debris) and with lentic local conditions.VES included dip net sweeps to assess early life history stages (egg, larvae).We listened for calls of breeding males while slowly driving or walking roads adjacent to wetlands, riparian areas and large water bodies or while floating the Yellowstone River.For road calling surveys, we stopped vehicles every 0.8 -1.6 km in low-lying areas suitable for containing pooled water.We listened for bullfrog calls for 5 minutes at each stop.At accessible stops where bullfrog calls were heard, we used VES to confirm adult presence and looked for egg masses and larvae.At inaccessible sites, we only detected adult males aurally.
2012 and 2013.To permit inference about bullfrog occupancy in the Yellowstone River floodplain, we used a generalized random tessellation stratified design to select survey sites.Strata were lacustrine, palustrine, and riverine wetland types (Cowardin et al. 1979) within the Yellowstone River floodplain.These strata describe >98% of aquatic habitats present outside of the main channel.We defined a site as an isolated lentic habitat (e.g., pond) or lotic environment with lentic local conditions (e.g., backwaters and sidechannels).This approach resulted in 102 sites that we were able to access in 2012.We surveyed 97 of these 102 sites in 2013.The up-and down-stream extent of surveys were identical in 2012 and 2013.Thirty sites that we surveyed in 2010 were also surveyed in 2012 and 2013.The downstream extent of 2010 surveys was ≈ 60 km greater than 2012 and 2013 surveys.
In 2012, we surveyed sites 1 -4 times between June 25 and August 24.We allowed at least 14 days between repeat surveys of a site.In 2013, we surveyed our whole pool of sites twice: once between July 10 -19 and once between August 10 -28.We surveyed all habitats between the hours of 8:00 and 20:00 using VES and calling surveys.In areas where we detected postmetamorphic bullfrogs, we used VES and dip net surveys to determine presence of egg and larval bullfrogs.Only presence/not-detected data were used for analyses.
During each survey, we quantified variables that may affect bullfrog detection: year, date, time of day, weather, water temperature, dissolved oxygen and maximum depth (Fuller et al. 2011;Ultsch et al. 1999).Weather was noted as clear, partly cloudy, overcast or rain.Water temperature and dissolved oxygen were measured once per site approximately 50 cm from shore at mid-depth in the water column with a YSI Professional Plus Multiparameter Meter (Yellow Springs, OH).We measured water depth to the nearest 0.1 m with a stadia rod at the deepest location we found in each site.
We used ArcGIS (v.10) with a National Wetlands Inventory layer (U.S. Fish and Wildlife Service 2010) to quantify covariates predicted to affect bullfrog occupancy: distance from public access sites, distance from human-modified habitats, habitat area, water depth, percent cover of emergent vegetation and wetland habitat type.We provide rationale for these covariates in the Data Analysis section.
We estimated the shortest Euclidean distances between the nearest edges of sites.We defined public access sites as any federal, state, county or city fishing access site, boat ramp and riverside public park.We identified these sites while conducting surveys and by using road maps and state, county and city websites.To identify human-modified habitats, we used our survey observations and wetland habitat codes [spoil] and x [excavated]) from the National Wetlands Inventory layer (Cowardin et al. 1979).
To estimate habitat area, we walked the perimeter of the surveyed habitat with a GPS (Montana 650t v., Garmin Ltd. 1996-2013) and calculated the area of the resulting polygon.To characterize water permanency, we described water depth at the deepest measured point and used the shallowest value for sites that were surveyed multiple times.We used visual surveys to estimate the percent cover of emergent vegetation.

Data analysis
To describe bullfrog spread, we compared the observed up-and downstream extent of bullfrog breeding and any bullfrog life stage in 2010 to the distal detections in 2012 and 2013.We built models to examine how bullfrog occupancy (ψ) in 2012 and colonization (γ) in 2013 were related to 3 groups of variables: habitats we suspected were more likely to have human-mediated introductions due to high use and access, habitats that were anthropogenic or where conditions had been substantially modified, and habitats with permanent waters.We separately assessed model support for presence of bullfrog breeding and any bullfrog life stage.
We posited that ψ and γ should increase with proximity to sites with high potential for people to introduce bullfrogs.Models describing potential human-mediated introductions include the covariates: access (public-access site or not) and distance of a surveyed site to the nearest publicaccess site (distance to public-access site).We also posited that ψ and γ should be high where humans have created habitat conditions suitable to bullfrogs such as deep, lentic sites with permanent hydroperiods (Fuller et al. 2011).Models describing human-modified habitats include the covariates: habitat type (human-modified habitat or not) and distance of a surveyed site to a humanmodified habitat (distance to human-modified habitat).Finally, we posited that ψ and γ should be high in stable, permanent waters that provide overwintering habitat and are less vulnerable to scouring flows or late summer desiccation (Fuller et al. 2011).Models describing stable, permanent habitats include the covariates: habitat area, water depth and percent cover of emergent vegetation (emergent).
We also assessed support for models that combined covariates of human-mediated introductions or human-modified habitats with covariates of stable, permanent water, for a null model with only the intercept and for a global model with all site-specific covariates.Prior to testing models, we transformed continuous covariates to Z scores with a mean of 0 and standard deviation of 1 (Donovan and Hines 2007).We did not find that covariates were correlated with one another (r ≥ 0.7).
To compare models, we used the original parameterization of the single species, multiseason model in PRESENCE 4.0 (MacKenzie et al. 2003).This approach allows for missing data and estimates site ψ, γ, extinction (ϵ) and detection (p) probabilities as functions of site and survey-specific covariates.We began by fitting all of the hypothesized models for p while using a global model structure for ψ and γ that included all site-specific covariates.We modeled ϵ as constant because we failed to detect bullfrogs in only 1 site in 2013 where we detected bullfrogs in 2012, which suggests ϵ is low.We compared detection models with combinations of our covariates of detectability.We ranked models with Akaike's Information Criterion modified for small sample size (AIC c ) and used the best ranked model (lowest AIC c ) to test occupancy and colonization models.We evaluated the strength of evidence for each of our hypothesized models based on AIC c and the resulting Akaike weights (w i ) for each model (Burnham and Anderson 2002).We used a ∆AICc < 4 as adequate separation to indicate supported models, but still considered models with a ∆AICc < 7 as having some support (Burnham et al. 2011).We used model averaging among the top models (AIC c < 7) to test if covariate coefficients differed from 0.

Results
We detected bullfrog breeding in 10 sites in 2010 and 12 sites in 2012 and in 45 sites in 2013 (Table 1, Figure 1).We detected any bullfrog life stage in 32 sites in 2010, 39 sites in 2012 and 58 sites in 2013 (Table 1, Figure 2).Of the 30 sites that we surveyed in all 3 years, we detected any bullfrog life stage in 1 site in 2012 and 7 sites in 2013 where we did not detect them in 2010.We never detected bullfrogs in 7 of the 30 sites and we always detected bullfrogs in 15 of the 30 sites across all years.
In 2013, we detected breeding bullfrogs 12 km upstream and 39 km downstream from where we detected them in 2010 (Table 1 and Figure 1).Changes in the up-and downstream extents of bullfrog breeding detections were greater between 2012 and 2013 surveys than between 2010 and 2012 surveys (Table 1).In 2013, we detected any bullfrog life stage 12 km upstream and 34 km downstream from where we detected them in 2010 (Table 1 and Figure 2).Changes in the up-and downstream extents of any bullfrog life stage detections were similar between 2012 -2013 surveys and between 2010 -2012 surveys (Table 1).
In 2012, we surveyed 102 sites and visited 60 of the 102 sites four times, 20 sites three times, 10 sites twice and 12 sites once.Naïve (unadjusted) detection of breeding and any life stage of bullfrogs were 100% and 91%, respectively.In 2013, we surveyed 97 sites twice.Naïve detection rates for breeding and any life stage of bullfrogs were 65% and 57%, respectively.

Breeding
No single model was clearly supported over others we considered (Table 2).Covariates that describe habitats with deeper waters and emergent vegetation were most supported.The coefficient of water depth was different from zero, with Ѱ increasing with water depth (Table 3, Figure 3).No covariate coefficients of γ differed from zero.We found similar support (∆AIC c < 4) for models that included the same predictor variables for Ѱ and γ and for models that only included predictor variables for Ѱ.The model that held Ѱ and γ constant had less support but was still plausible (∆AIC c < 6).
Models that held p constant in 2012 but varied p by survey in 2013 were most supported (∆AIC c > 8) and had a combined model weight of 95%.Probability of detection (± 95% C.I.) in the top model was 0.94 ± 0.04 for 2012 surveys, 0.72 ± 0.18 for July 2013 surveys and 0.26 ± 0.12 for August 2013 surveys.

Any bullfrog life stage
The model with the strongest support was Ѱ(public-access site + water depth) γ(.), indicating that occupancy was related to human-mediated introductions and stable, permanent waters (Table 2).The coefficient of public-access sites also occurred in the only other two models with ∆AIC c < 7. Coefficients of public-access sites and water depth were positive and different from 0 (Table 3).Ѱ was higher in sites with public access and in deeper waters (Figure 3).We found minimal support that γ was related to publicaccess sites (Table 2), but this coefficient did not differ from 0 (Table 3).Models that held detection probability constant in 2012 but varied it by survey in 2013 were most supported (∆AIC c > 9) and had a combined model weight of 98%.Probability of detection (± 95% C.I.) in the top model was 0.81 ± 0.09 for 2012 surveys, 0.50 ± 0.12 for July 2013 surveys and 0.83 ± 0.12 for August 2013 surveys.

Discussion
Our survey data describe an invader that is rapidly spreading along the Yellowstone River floodplain.As of 2013, bullfrogs occupy at least 58 sites along 107 km of floodplain.The number of sites with evidence of bullfrog breeding increased from 12 in 2012 to 45 in 2013.These breeding data represent a spread of the bullfrog population rather than an artefact of increased sampling effort because the number of surveyed sites was similar across years (Table 1), sampling intensity was greater in 2012 than in 2013, and we detected bullfrogs in 2012 and 2013 in 7 of 30 sites where we failed to detect them in 2010.These data also indicate that bullfrogs are firmly established in the Yellowstone River floodplain and seem poised to rapidly spread to uninvaded habitats.
Our observations of spread indicate that bullfrogs are dispersing in up-and downstream directions (Figure 2).However, most occupied sites in 2012 and colonized sites in 2013 were within the distribution range observed in 2010.This pattern may indicate that successful colonization events in our area tend to be associated with shorter distances.While not included in our surveys, it is also possible that bullfrogs from sites in adjacent uplands may also be contributing to the spread within the Yellowstone River floodplain.
We found that bullfrog spread was greater in a downstream direction.Since 2010, bullfrogs have spread 34 km downstream and only 12 km upstream.In addition to active movement along the floodplain, it is possible that river flow may facilitate downstream bullfrog spread.Indeed, flow can facilitate and accelerate the downstream spread of multiple invasive species, including dreissenid mussels (e.g., Johnson and Padilla 1996) and crayfish (Bubb et al. 2004).We observed bullfrogs along the main channel and in active side channels, so there is opportunity to disperse passively during high flows or actively downstream during lower flows.In fact, larvae in lab conditions respond to current by moving downstream (Schmidt et al. 2011).
Another potential driver of bullfrog spread is human-mediated introductions.Humans have introduced bullfrogs for food, recreational hunting, bait and pest control and as released pets no longer suitable for aquaria (Bury and Whelan 1984;Jennings and Hayes 1985).Given the frequency that invasive species are intentionally and accidently transferred to new waters (e.g., Johnson et al. 2001) and the positive association between bullfrog Ѱ and public-access sites, it is possible that bullfrogs have been introduced to multiple public-access sites in the Yellowstone River region.Publicaccess sites are especially vulnerable to aquatic invasive species introductions because they receive high traffic and are easily accessible (Johnson et al. 2001).Therefore, efforts to prevent further introductions should include education and outreach campaigns that target public-access users.Public-access sites may also be key areas to focus control efforts.
We found that bullfrog eggs and larvae were associated with deeper waters that have emergent vegetation, while adults were associated with deeper waters (Table 2 and Figure 3).Deep waters characterize habitats with permanent water, which are essential in this region because larvae must overwinter before reaching metamorphosis (Adams and Pearl 2007;Cook and Jennings 2007;Peterson et al. 2013).In colder climates, like Montana, adults may also require deeper waters for overwintering because they cannot tolerate freezing conditions.Finally, deeper waters may provide important refuge from high river flows that can scour shallower habitats (Fuller et al. 2011).For these reasons, draining permanent waters has been recommended and successfully applied by others to reduce bullfrog densities (e.g., Doubledee et al. 2003;Maret et al. 2006;Adams and Pearl 2007;Fuller et al. 2011).
Indeed, 2013 appears to be a period of bullfrog expansion in our study region -we observed a large increase in the number of breeding sites relative to 2012.This increase may have been caused by antecedent flow conditions, or lag effects in bullfrog abundance or breeding detectability.High flows are predicted to disturb habitat conditions, like emergent vegetation, that favor bullfrogs in active floodplains (Fuller et al. 2011).In 2011, the Yellowstone River had a 29year flow event (2,006 m 3 •s -1 ) that activated the floodplain, including many sites that we surveyed.Of the sites surveyed in 2010 and 2012, bullfrogs were detected both years at all 4 sites that were isolated from the main channel, but only at 4 of 9 sites that occurred within the floodplain.It is likely the bullfrog distribution and abundance in 2012 were influenced by the 2011 high flows.Additional hypotheses are that there was an inherent lag effect in bullfrog abundance and dispersal that was overcome in 2013 (sensu Crooks and Soulé 1999) or more optimal breeding conditions during the 2012 summer (compared to 2011) increased breeding detectability in July 2013.
We found that sampling occasion influenced detection of bullfrog breeding and any bullfrog life stage, similar to other studies (e.g., Gooch et al. 2006).In 2013, detection of bullfrog breeding was high in July but low in August while the opposite was observed for juvenile and adult stages.Sites where we detected breeding in July 2013 may have had poor egg and larvae survival.Differences in p could also be related to reduced hydroperiods in 2013; August discharge levels in 2013 (48 m 3 •s -1 ) were lower than 2012 (67 m 3 •s -1 ; USGS station 06214500).In response to less surface water in 2013, larvae may have metamorphosed between sampling occasions or moved to deeper waters in habitats that were connected to the main channel.However, the low flows in August 2013 appeared to concentrate juvenile and adult bullfrogs, thereby increasing their detectability.

Conservation and management
Our study identified public-access sites and permanent waters as habitats on which to focus bullfrog prevention, early detection, direct removal and restoration.Nevertheless, large rivers like the Yellowstone are difficult places to control introduced species because of their size, complex habitat and access (Sepulveda et al. 2013;Tyus and Saunders III 2000).In regulated rivers, one strategy for bolstering native assemblages and suppressing introduced species is to restore the natural flow regime (e.g., Marchetti and Moyle 2001).In undammed rivers like the Yellowstone, maintaining natural flow variability is important, but direct removal of introduced species from native strongholds may also be necessary (Fuller et al. 2011;Propst et al. 2008;Tyus and Saunders III 2000).For bullfrogs, the effectiveness of direct removal efforts at the scale of the Yellowstone invasion is questionable because bullfrogs are evasive, have strong density-dependence in both the larval and post-metamorphic stages, and can quickly recolonize from neighbouring wetlands (Adams and Pearl 2007;Doubledee et al. 2003;Govindarajulu et al. 2005).

Figure 1 .
Figure 1.Survey sites where bullfrog breeding was found (dark circles) and not detected (gray circles) in 2010, 2012 and 2013 in the Yellowstone River region near Billings, MT.The Yellowstone River flows southwest to northeast.Table S1 provides geo-referenced records of bullfrog breeding detections.

Figure 3 .
Figure 3. Association of mean site-specific parameter estimates (± 1 SE) of occupancy by bullfrog breeding life stages (gray bars) and any life stage (dark bars) with water depth (m) at the deepest location found in each site.

Table 1 .
The number of sites surveyed (N), the number of sites where we detected bullfrog breeding and any life stage and the upstream and downstream extent (river km) of bullfrog breeding and any life stage in the Yellowstone River floodplain in 2010, 2012 and 2013.

Table 2 .
Best supported models predicting bullfrog breeding and any life stage occupancy (ψ) and colonization (γ) in the Yellowstone River floodplain in 2012 and 2013.The model structure for detection probability was p (2012, July 2013, August 2013) and we held this model structure constant across all ψ models.All models also included the intercept.K is the number of parameters, AICc is the second-order Akaike information criterion, AICc Wt is the second-order Akaike weight, and deviance is -2(loglikehood) of each model.Only models with ΔAICc < 7, the global model and the null model are shown.

Table 3 .
Model -average coefficients ( ), 95% confidence intervals (CI), and cumulative Akaike weights (AICc Wt) for variables predicting bullfrog breeding and any life stage occupancy (ψ) in 2012 and 2013.Coefficients for detection are from the top-ranked model.Only variables included in models with ΔAICc < 7 are shown.Dashed lines (-) indicate the variable was not including in models with ΔAICc < 7.