Distribution, Spread, and Habitat Predictability of a Small, Invasive, Piscivorous Fish in an Important Estuarine Fish Nursery

Invasive species often cause negative ecological and economic effects. Florida has >20 established invasive fish species but only seven exist in saltwater. The present study examined Belonesox belizanus (Pike Killifish), a Central American euryhaline fish introduced to Tampa Bay (west-central Florida) in the early 1990s, which has quantifiably reduced populations of small-bodied native fishes and may compete with prized sportfish juveniles in estuarine nursery habitat. Long-term monitoring revealed that B. belizanus occurs in estuarine waterbodies along a 31-km stretch of the bay’s eastern fringe, with a second, smaller population in two western tributaries. Spread rate was estimated to be 5.5–13 km year−1, intermediate among invasive poeciliids. A novel implementation of boosted regression tree modeling to assess B. belizanus habitat predictability found greater probability of presence with decreasing water depth and pH, whereas presence tended to be greatest at polyhaline salinity. It is hypothesized that B. belizanus distribution in Tampa Bay is constrained by deep, seawall habitats acting as ecological barriers. Further B. belizanus spread therefore may be most likely to occur by human release (from aquaria or bait buckets) or bird carry-off. Newly restored tidal habitat within the current range probably will be invaded quite quickly by B. belizanus.


Introduction
Invasive species, habitat degradation or loss, and pollution are the most important threats to native fishes in the USA [1].Invasive species are defined by the United States government as "with regard to a particular ecosystem, a non-native organism whose introduction causes or is likely to cause economic or environmental harm, or harm to human, animal, or plant health" [2].The number of established nonnative fishes in the USA more than doubled between 1970 and 2013 [3].Negative ecological effects from invasive fishes can include predation, competition, and interference [4].Economic losses from invasive fishes have been estimated at several billion US dollars annually [5].Successful establishment of invasive fishes in new habitats, whether through intentional or unintentional releases, may depend on a variety of factors including broad physiological tolerance, rapid dispersal, and similarity of habitat to the native range [6].
Prediction of invasive fish presence and habitat use is one of the most important applications of species distribution modeling [7].Establishing species-habitat relationships that are general and broadly applicable rather than idiosyncratic to a narrow set of conditions is desirable for allowing successful prediction of presence in unsampled areas [8].Prediction serves as an important means of validating the generality of observed species-habitat relationships and was recently argued by Houlahan et al. [9] to be essential for demonstrating ecological understanding.Peninsular Florida is likely to contain more invasive freshwater fish than any other similarly sized area in the world [10], probably because the warm climate facilitates survival of tropical fish that are released from freshwater aquaria or that escape from fish farms [11].Although freshwater invasive fishes are prevalent, there are relatively few established invasive saltwater or estuarine fishes in Florida: Schofield and Loftus [12] described 22 permanently established freshwater fishes, whereas Schofield et al. [13] noted only two permanently established marine fishes (the lionfishes Pterois volitans and Pterois miles).Six of the permanently established freshwater fishes are euryhaline and occur in saltwater habitat [13].The Pike Killifish Belonesox belizanus Kner (Pisces: Poeciliidae) is among these species and was introduced to Miami-Dade County, south Florida, in 1957 when ~50 excess medical research specimens were released into an agricultural canal [14].Although the main invasive distribution of B. belizanus is in south Florida (Figure 1), a secondary introduction occurred in the Tampa Bay watershed (west central Florida), possibly as a result of escape from a tropical aquaculture facility [15], and is the subject of the present study.
B. belizanus is native to Central American Atlantic drainages south of the Punta del Morro (Mexico) to northeastern Costa Rica [16,17] and although the largest species within the Poeciliidae, it is still relatively small (maximum size of 200-mm total length; [17]).All life stages of B. belizanus are carnivorous [17], with prominent teeth and a novel gape adaptation that allows the species to consume prey that are relatively large for its size [18].The primary prey for all life stages (including neonates) is fish [14,[19][20][21].The ability to switch to invertebrate prey such as small shrimps is considered to allow B. belizanus to be a successful invader despite its specialized, invariant feeding mechanism [22].In south Florida, B. belizanus was observed to deplete small resident fish populations after introduction [14,19] and the species was recently classified at the upper end of medium risk for potential negative effects [23].Recognition of the ecological risk from B. belizanus is evidenced by the species being prohibited for possession in a number of American states, although not Florida.
There is strong quantitative evidence that the Tampa Bay population of B. belizanus has significantly reduced abundance of small-bodied resident fishes, in particular Sailfin Molly Poecilia latipinna, Eastern Mosquitofish Gambusia holbrooki, and Sheepshead Minnow Cyprinodon variegatus [24].Although there is little evidence for direct effects on abundance or size distribution of estuary-dependent migratory species, which generally grow to larger sizes and recruit from larger source populations than resident fishes, the depletion of resident prey fishes by B. belizanus is of concern for its potential indirect effects on estuary-dependent species through competition [15].In particular, Tampa Bay provides important nursery habitat for Common Snook Centropomus undecimalis, a prized sportfish supporting a substantial recreational fishery (>277,000 directed trips year −1 worth $15.8 M year −1 ; [25]), which experienced a considerable cold-kill-related decline in recent years [26] and occurs in estuarine areas in which B. belizanus is also found [27][28][29].
Knowledge of B. belizanus spread, distribution, and habitat associations in Tampa Bay is therefore of considerable importance both from the perspective of biodiversity (structuring of resident native fish assemblages) and potential effects on managed fisheries.Qualitative accounts of B. belizanus habitat in the native range describe occupation of a variety of shallow-water, low-flow habitats ranging from freshwater to hypersaline, including mangrove swamps, sinkholes, ponds, creeks, and river backwaters [17,[30][31][32].Quantitative evaluations of B. belizanus habitat have included multispecies analyses in the native range [33][34][35][36] and a detailed single-species evaluation of associations with physicochemical parameters from several datasets collected in south Florida [37].The latter study found B. belizanus density was correlated with salinity and temperature at the southernmost, estuarine sites, whereas pH was most associated with changes in density at the northern, inland sites.The present study used long-term monitoring and shorter duration sampling efforts [24] to ask three main questions.over a greater spatial area of Florida and other states in the future [38].Third, what are the habitat associations of B. belizanus?Following the philosophy of Houlahan et al. [9] that prediction should be used to aid ecological understanding, the present study employed a novel method to assess predictability of habitat use in test waterbodies based on boosted regression tree models developed from independent, training waterbodies.This allowed examination of the extent to which observed B. belizanus-habitat relationships were general and broadly applicable rather than idiosyncratic and limited to specific waterbodies.The present study is the first quantitative, single-species examination of B. belizanus-habitat relationships based solely on data collected in situ at the time of sampling, and demonstrates a method that can also be used for other species and locations.The main conclusions of the present study are that the principal Tampa Bay population of B. belizanus is distributed in eastern estuarine waterbodies, with a smaller, spatially restricted western population that is likely to have been a separate introduction.Although the rate of B. belizanus spread is comparable to that of some other invasive fishes, the distribution may be constrained by soft ecological barriers, which emphasizes the need to manage human-mediated spread.Success of B. belizanus in Tampa Bay appears facilitated by its spread rate, broad environmental tolerance, and relatively distinct ecological niche.

Materials and Methods
All data used in the present study are provided as Supplementary Material.

Study Area and Sampling Techniques
Tampa Bay is Florida's largest open-water estuary (1031 km 2 ; [39]), a subtropical system that, although greatly altered by urbanization, still has relatively large areas fringed by mangroves, salt marshes, and other natural vegetation [29,40,41].Four main tributaries (the Hillsborough, Alafia, Little Manatee, and Manatee rivers) provide freshwater inflow, with approximately 100 smaller tributaries also draining to the bay (Figure 1; [42]).The smaller tributaries and creeks range from relatively natural mangrove-or saltmarsh-lined waterbodies to urbanized systems with considerable amounts of seawall and riprap.Physicochemical conditions vary widely and salinity ranges from nearly fresh to polyhaline, depending on tides, precipitation, and wind.Tampa Bay is described in detail by Lewis and Estevez [39].
Analyses presented herein used data collected in 1989-2008 during diurnal shoreline seine sampling within estuarine waterbodies of the Tampa Bay watershed by the Florida Fish and Wildlife Conservation Commission, Fish and Wildlife Research Institute's Fisheries-Independent Monitoring program [24].Data from trawling in deeper areas and seines deployed away from the shore were not considered because B. belizanus was never collected in these habitats.In all sampling, B. belizanus specimens were identified, counted, and measured (standard length, mm).Descriptive habitat data and physicochemical parameters were also recorded, as described further below.Secchi depth as a measure of water clarity was also measured, but is not included in analyses because sampling sites were sufficiently shallow that the Secchi disk was almost always visible on the substrate.Further details of the sampling procedures can be found elsewhere [24,43].

Distribution and Spread Rate
Coordinates of sites at which B. belizanus was present or absent during seine sampling were overlaid on maps of the Tampa Bay watershed to assess the species' distribution.Long-term sampling in the lower Little Manatee River watershed was used to examine the rate of spread following first occurrence.This watershed was selected to examine rate of spread because it offered the most comprehensive, consistent longitudinal sampling and monthly sampling had already commenced before first occurrence (in contrast to the Alafia River, for example, wherein sampling was only seasonal at the time of first occurrence).Monthly 21.3-m shoreline seine sampling in the Little Manatee River watershed occurs in the lowermost 13 km of the river main stem and accessible adjacent habitats (Figure 2).Sampling is stratified into upstream and downstream zones 7 km upstream of the river mouth (Figure 2), with 4 samples per month in each zone.The temporal sequence of observations suggested that B. belizanus invaded the Little Manatee River from downstream to upstream.Following Ben Rais Lasram et al. [44], the rate of spread (S) was calculated: where d is distance upstream of the mouth of the Little Manatee River (km), D is number of days from the first observation in the river, t is a given observation of B. belizanus, and t − 1 is the previous downstream observation of B. belizanus.
Fishes 2017, 2, 6 5 of 23 occurrence.This watershed was selected to examine rate of spread because it offered the most comprehensive, consistent longitudinal sampling and monthly sampling had already commenced before first occurrence (in contrast to the Alafia River, for example, wherein sampling was only seasonal at the time of first occurrence).Monthly 21.3-m shoreline seine sampling in the Little Manatee River watershed occurs in the lowermost 13 km of the river main stem and accessible adjacent habitats (Figure 2).Sampling is stratified into upstream and downstream zones 7 km upstream of the river mouth (Figure 2), with 4 samples per month in each zone.The temporal sequence of observations suggested that B. belizanus invaded the Little Manatee River from downstream to upstream.Following Ben Rais Lasram et al. [44], the rate of spread (S) was calculated: where d is distance upstream of the mouth of the Little Manatee River (km), D is number of days from the first observation in the river, t is a given observation of B. belizanus, and t − 1 is the previous downstream observation of B. belizanus.

Habitat Predictability
Habitat predictability of B. belizanus was investigated using data collected monthly in 2006 or 2007 with 9.1-m seines in 11 estuarine waterbodies within the species' range in the Tampa Bay watershed, consisting of small creeks, canals, and adjacent non-creek/canal habitats (river main-stems or bayous) (Figure 1, Tables 1 and 2; [24]).The original purpose of sampling in these waterbodies was to assess the value of tidal creeks/canals as nekton habitat, with the adjacent non-creek/canal habitats generally also sampled to provide a contrast.Each waterbody was divided into longitudinal spatial strata (~200 m long) and daylight sampling was undertaken monthly at randomly chosen locations between each creek/canal's mouth and the upstream extent that could be reached within time or access constraints.Sites were chosen from a set of predefined geographic coordinates with a random number generator.The number of samples per waterbody was proportional to the extent of the waterbody that was sampled, with two samples per spatial stratum per month; thus, the number of seine samples ranged from four per month in smaller waterbodies (250-300 m long) to 10 per month in larger waterbodies (1000 m long).Adjacent waterbodies (rivers or bayou) also received two samples per month per stratum, with one stratum located upstream and downstream of the creek/canal mouths.A total of 935 shoreline seine samples was collected.Additional detailed description of the studied waterbodies is provided in [21,45].
Relationships between B. belizanus probability of presence and a variety of habitat predictor variables were examined using Boosted Regression Trees (BRT; [46]).As described by Elith et al. [46], BRT is a machine learning method that combines regression trees (models that relate a response to their predictors by recursive binary splits) and boosting (an adaptive method for combining many simple models to give improved predictive performance), resulting in final models that are essentially additive regression models in which individual terms are simple trees, fitted in a forward, stagewise fashion.The advantages of BRT include no need for prior data transformation or elimination of outliers, ability to fit complex nonlinear relationships, and inclusion of interactions between predictors [46].
Habitat variables were assessed in the BRT modeling for their potential to influence the probability of presence of B. belizanus based on previously observed or hypothesized relationships.Julian day was included to represent seasonality [37].Waterbody type was included to test for the potential importance of mesohabitat, given previous qualitative observations of habitat preference (Miller [30]: "swamps and the eddies of streams"; Turner and Snelson [20]: " . . .habitats with some degree of cover, such as mangrove swamps and narrow, shallow canals with dense aquatic vegetation . . .not found in large, open canals or borrow ponds"; Greenfield and Thomerson [31]: " . . .shallow water and areas of heavy vegetation", "temporary ponds that form during . . .flooding", "larger rivers and streams . . .usually found close to shore, in slower-moving backwater areas"; and Miller et al. [17]: " . . .ciénegas, lagoons, ponds, and slow-moving parts of creeks and rivers . . ."). Waterbodies were classified as creek (natural tributaries of larger rivers), canal (dredged inlets of larger rivers), bayou (large open embayment of a river), and river (mainstems of the two main rivers within the B. belizanus range, Alafia and Little Manatee).Dominant shoreline type at each sampling site was included because of qualitative observations of differences in B. belizanus presence between vegetation type [20], and was characterized as emergent marsh (primarily black rush, Juncus roemarianus), mangrove (red mangrove, Rhizophora mangle, or black mangrove, Avicennia germinans), small terrestrial vegetation (primarily leather fern, Acrostichum danaefolium), shrubs/trees, or artificially hardened (e.g., sea wall or rip rap).Associated with shoreline type were qualitative variables for whether or not the shoreline was inundated during sampling and whether or not the shoreline vegetation was overhanging the sampling site; these variables could influence the probability of presence along shoreline habitat, as shown for other fishes [47,48].Substrate type was represented in the BRT models as mud, mud-sand, sand, or some oyster/rock.Water depth (at the seine's bag, i.e., 1.5 m from shore; with sampling limited to sites 1.2 m deep or less) and site gradient (calculated as the difference in water depth between the seine's bag and the seine's wing at the shore, divided by the distance between bag and shore, and expressed as a percentage) were included to address these variables' potential importance for predation risk [49].Physicochemical variables previously shown to correlate with B. belizanus density in south Florida [37] were also included in the BRT models: salinity, temperature, dissolved oxygen, and pH.These physicochemical variables were measured at 0.2-m depth with a Hydrolab ® multiprobe (OTT Hydromet, Kempten, Germany), reflecting B. belizanus's tendency to inhabit surface waters [30].In order to examine habitat predictability for B. belizanus presence, a novel approach based on creating two sets of BRT models was implemented.Hereinafter the models are referred to as the fixed subset and randomly selected subset models.For both types of models, training data are defined as the data used to formulate the model, whereas test data were not used for model formulation and instead were used to test the models' predictive performance, i.e., the test data were independent of the training data.For the fixed subset model, test data were represented by the data collected in three waterbodies (Andrews Creek, Curiosity Creek, and Question Mark Creek; n = 240 samples, 25.7% of the 935 total samples), which were chosen on the basis of their habitat characteristics being within the range encompassed by the other eight waterbodies (Tables 1 and 2).Training data thus consisted of the 695 samples (74.3%) from these eight other waterbodies.The fixed subset BRT model was implemented in the dismo package [50] of the R software environment [51] using a bag fraction (i.e., the proportion of data to be selected at each step of the model) of 0.5, a learning rate (i.e., the shrinkage parameter determining the contribution of each tree to the growing model) of 0.005, and tree complexity of 2 (i.e., allowing first-order interactions) in order to achieve Elith and coworkers' [46] recommendation of >1000 trees for reliable model predictions.Exploration of a higher order of tree complexity (5) did not yield substantial improvement in the model's predictive performance, so tree complexity of 2 was implemented for the final analysis.Given the stochasticity in BRT model formulation [46], 1000 iterations of the fixed subset model were undertaken in order to examine the variability in model predictive performance and the consistency in relationships between B. belizanus presence and habitat.Model fitting, tenfold crossvalidation, and prediction of the model to the test data were undertaken for each of the 1000 iterations of the fixed subset model.
The randomly selected subset BRT model consisted of randomly selecting (without replacement) 695 samples of the total 935 samples to use as training data, then repeating the steps described above for the fixed subset model: model fitting, crossvalidation, and prediction to the test data (i.e., the randomly remaining 240 samples).This whole procedure was repeated 1000 times, again to illustrate the variability in model predictive performance and consistency in B. belizanus presence-habitat relationships.
The predictive performance of the fixed and randomly selected subset BRT models in correctly modeling presence of B. belizanus was assessed using the area under the receiver operating characteristic curve (AUC).As noted by Huang and Frimpong [7], AUC values are commonly interpreted as follows: AUC > 0.9 is excellent; 0.8 < AUC < 0.9 is good; 0.7 < AUC < 0.8 is fair; 0.6 < AUC < 0.7 is poor; and 0.5 < AUC < 0.6 is failed (i.e., similar to random guessing).Huang and Frimpong suggested that this classification scheme is too demanding for prediction to test data, given that the scheme was developed for the purpose of assessing model fit rather than transferability (prediction); instead, they used a threshold of AUC < 0.6 to indicate limited transferability from training to test data for prediction [7].This threshold was also considered appropriate for the present study.Similar to the approach used by Broennimann and Guisan [52], the variability in AUC over the 1000 iterations of the fixed subset and randomly selected BRT models allowed 95% intervals to be calculated, based on the 2.5th and 97.5th percentiles of AUC.The overlap of the 95% intervals was assessed to examine the evidence for differences in predictive performance between fixed subset and randomly selected BRT models for training, crossvalidation, and prediction to test data.
In addition to assessing predictive performance with AUC, the ability of the fixed subset BRT to predict presence in the three test waterbodies was assessed by summing the probability of presence over all samples for each waterbody to give an overall predicted % presence for each iteration of the model, with 95% intervals calculated based on the 2.5th and 97.5th percentiles of the 1000 iterations.The predicted estimates of % presence were compared to the actual % presence in each waterbody, for which 95% intervals were generated by bootstrapping-randomly sampling with replacement, for the same sample size as the original sample [53]-the presence/absence data 1000 times, then calculating the 2.5th and 97.5th percentiles of the % presence from each of the 1000 resamples.The relative contribution [54] of the habitat predictor variables to B. belizanus presence was compared between the fixed and randomly selected subset BRT models using the 95% intervals from the 1000 model iterations.The consistency of B. belizanus presence-habitat relationships between the fixed and randomly selected subset models was assessed with partial dependence plots from R's gbm package [55].The partial dependence plots illustrated the marginal effect of each habitat predictor on B. belizanus probability of presence after accounting for the average effects of all other predictors in the models, and were based on the 95% intervals of the marginal effects from the 1000 model iterations.

Distribution and Spread Rate
The first specimen of B. belizanus collected in the Tampa Bay watershed was made on 21 February 1994 from a canal off the Alafia River (at 3.2 km from the river mouth; Figure 1).During extensive sampling occurring in 2006/2007, B. belizanus was mostly collected in a number of waterbodies draining either directly to eastern Tampa Bay, to the Alafia River, or to the Little Manatee River (Figure 1).Data collected following this period (2008-2016), primarily as a continuation of long-term monitoring for status and trends in estuary-dependent sportfishes and nekton communities, found B. belizanus in similar areas within this range [56,57].The straight-line distance from the northernmost to southernmost collection locations is approximately 31 km.B. belizanus was also collected in two small urban creeks in the city of St. Petersburg, on the western fringe of Tampa Bay (Figure 1).
The first collection of B. belizanus in the Little Manatee River watershed was made on 13 October 1998, approximately 5 km from the mouth of the river, with the farthest upstream collection within the sampling area at approximately river km 13 of the main stem in October 2000 (Figure 2).Four collections between October 1998 and November 1999 gave estimated spread rates of 15-35 m day −1 (5.5-13 km year −1 ; Figure 2).

Habitat Predictability
B. belizanus was found in 7% of shoreline seine samples from the 11 estuarine waterbodies studied during 2006/2007, with the greatest frequency of presence (12%-14%) in two polyhaline mangrove-dominated creeks (Andrews Creek and Cockroach Creek), a mesohaline marsh-dominated creek (Wildcat Creek), and a mesohaline canal mostly lined with shrubs/trees (Question Mark Creek) (Table 3).The species occurred along all shoreline and bottom types sampled with seines, across a broad range of physiochemical variables: salinity 0.2-29.3,temperature 11.8-30.1 • C, dissolved oxygen 0.5-14.7 mg l −1 , pH 6.9-8.4,depth 0.1-1.0m, and gradient 0%-66.7%.The predictive performance of the BRT models for B. belizanus presence was greatest for the training data, with slightly more variation in the area under the receiver operating characteristic curve (AUC) for the 1000 iterations of the randomly selected data subsets (95% interval: 0.92-0.98)than for the fixed subsets (0.92-0.96) (Figure 3).The 95% interval of AUC for the crossvalidation of the training data overlapped for the randomly selected and fixed subsets of data, again with more variation for the randomly selected subset.With respect to prediction of test data that were not included in development of the models, there was significantly lower AUC for the fixed subset model compared to the randomly selected model, with the latter having broad variation which overlapped the AUC of the crossvalidation models (Figure 3).In general, the BRT models based on the randomly selected and fixed subsets of the data suggested a similar relative ranking of habitat predictor variable importance in predicting B. belizanus presence (Figure 5).The range of variability in importance across the 1000 iterations of the two BRT models was considerably greater for the model based on randomly selected subsets of the data.Julian day was the most important predictor, followed by dissolved oxygen.Depth was suggested to be as important as pH, salinity, temperature and shoreline type in the model based on randomly selected subsets of the data, whereas shoreline type in particular had greater importance in the model based on the fixed subset of data.Substrate, gradient, and waterbody had relatively little importance (generally 5% or less), and shoreline inundation and vegetation overhang had minimal importance (Figure 5).The BRT model for the fixed subset of the data provided a reasonable prediction of the presence of B. belizanus in one of the three test waterbodies: the 95% interval of the prediction for Curiosity Creek was 6.8%-7.4%,which was somewhat greater than the actual observed 5.0% presence, but well within the 95% interval for actual presence (Figure 4).In contrast, the BRT model underpredicted the presence in the two other waterbodies, and in the case of Andrews Creek, the 95% intervals for the predicted and actual presence did not overlap (Figure 4).In general, the BRT models based on the randomly selected and fixed subsets of the data suggested a similar relative ranking of habitat predictor variable importance in predicting B. belizanus presence (Figure 5).The range of variability in importance across the 1000 iterations of the two BRT models was considerably greater for the model based on randomly selected subsets of the data.Julian day was the most important predictor, followed by dissolved oxygen.Depth was suggested to be as important as pH, salinity, temperature and shoreline type in the model based on  In general, the BRT models based on the randomly selected and fixed subsets of the data suggested a similar relative ranking of habitat predictor variable importance in predicting B. belizanus presence (Figure 5).The range of variability in importance across the 1000 iterations of the two BRT models was considerably greater for the model based on randomly selected subsets of the data.Julian day was the most important predictor, followed by dissolved oxygen.Depth was suggested to be as important as pH, salinity, temperature and shoreline type in the model based on randomly selected subsets of the data, whereas shoreline type in particular had greater importance in the model based on the fixed subset of data.Substrate, gradient, and waterbody had relatively little importance (generally 5% or less), and shoreline inundation and vegetation overhang had minimal importance (Figure 5).The relationships between B. belizanus probability of presence and habitat predictors showed some similarities and some differences between the BRT models based on the randomly selected and fixed subsets of the data.Both the randomly selected and fixed subset BRT models suggested that the probability of presence increased from around Julian day 240 (i.e., September) until the end of the year, although the increase was less pronounced for the randomly selected model (Figure 6).A major difference between the randomly selected and fixed subset models was the representation of the effect of dissolved oxygen in hypoxic conditions (<2 mg l −1 ), for which the fixed subset model suggested a very high marginal effect on probability of presence (up to nearly 0.14) whereas the randomly selected subset model gave a marginal effect of around 0.04.This difference principally arose from high frequency of presence (0.6, i.e., 60% of samples) for B. belizanus during hypoxic conditions in Cockroach Creek in September 2007, with these observations being present in all iterations of the fixed subset model; the predicted effect of dissolved oxygen was more consistent between the randomly selected and fixed subset models over the rest of the variable's range.
Other habitat predictors with relatively high contributions for which the predicted marginal effect on B. belizanus presence was reasonably consistent between randomly selected and fixed subset BRT models were depth and pH, for which the probability of presence decreased by two to three times over the range of each variable (Figure 6).The predicted marginal effect of salinity was not mononotic and predicted B. belizanus presence as a peak between salinity of about 20 and 30 for both randomly selected and fixed subset models, whereas somewhat higher probability of presence at low salinity (<7) was only suggested by the fixed subset model.The increased probability of presence with increasing temperature, particularly above ~19 °C, was consistent for fixed and randomly selected BRT models.
The fixed subset model suggested some significant differences between shoreline types (with marsh having greater probability of B. belizanus presence than other types), substrate (greater The relationships between B. belizanus probability of presence and habitat predictors showed some similarities and some differences between the BRT models based on the randomly selected and fixed subsets of the data.Both the randomly selected and fixed subset BRT models suggested that the probability of presence increased from around Julian day 240 (i.e., September) until the end of the year, although the increase was less pronounced for the randomly selected model (Figure 6).A major difference between the randomly selected and fixed subset models was the representation of the effect of dissolved oxygen in hypoxic conditions (<2 mg L −1 ), for which the fixed subset model suggested a very high marginal effect on probability of presence (up to nearly 0.14) whereas the randomly selected subset model gave a marginal effect of around 0.04.This difference principally arose from high frequency of presence (0.6, i.e., 60% of samples) for B. belizanus during hypoxic conditions in Cockroach Creek in September 2007, with these observations being present in all iterations of the fixed subset model; the predicted effect of dissolved oxygen was more consistent between the randomly selected and fixed subset models over the rest of the variable's range.Other habitat predictors with relatively high contributions for which the predicted marginal effect on B. belizanus presence was reasonably consistent between randomly selected and fixed subset BRT models were depth and pH, for which the probability of presence decreased by two to three times over the range of each variable (Figure 6).The predicted marginal effect of salinity was not mononotic and predicted B. belizanus presence as a peak between salinity of about 20 and 30 for both randomly selected and fixed subset models, whereas somewhat higher probability of presence at low salinity (<7) was only suggested by the fixed subset model.The increased probability of presence with increasing temperature, particularly above ~19 • C, was consistent for fixed and randomly selected BRT models.
The fixed subset model suggested some significant differences between shoreline types (with marsh having greater probability of B. belizanus presence than other types), substrate (greater probability of presence over mud than oyster/rock or sand), and waterbody (greater probability of presence in canals and creeks than rivers, with bayou intermediate); however, the random subset model suggested a greater extent of overlap in the different levels of these predictors.A somewhat increasing probability of presence with increasing gradient was more pronounced in the random subset model than the fixed subset model (Figure 6).

Distribution and Spread Rate
The present study demonstrated that B. belizanus is well distributed in some estuarine portions of the Tampa Bay watershed, principally in the Alafia and Little Manatee River basins (Figure 1).From the first recorded presence in the Alafia River drainage in 1994, the species was subsequently detected in the Little Manatee River drainage in 1998.The mouths of these two rivers are approximately a straight-line distance of 17 km apart.Given the indentations along the Tampa Bay shoreline between the two rivers (Figure 1), as well as numerous convoluted potential inshore pathways (e.g., mosquito ditches and county road drainage ditches), the first presence of B. belizanus in the Little Manatee River nearly five years after the first recorded presence in the Alafia River is consistent with the estimated rate of spread in the Little Manatee River (5.5-13 km year −1 ; Figure 2).Subsequent sampling in non-tidal freshwater ditches that are ultimately connected to estuarine tributaries has shown that B. belizanus has spread extensively inland in these drainages [56].Collections of B. belizanus in south Florida following initial introduction indicates a spread rate of about 1.8 km in just over two months (from latitude 25.559516 • , longitude −80.331945 • in early December 1957 to approximately latitude 25.543571 • , longitude −80.331424 • on 18 February 1958; [14]).Extrapolation of this rate to a full year gives a value (9 km year −1 ) within the range estimated in the present study.Additional study of B. belizanus spread rates in south Florida based on the extensive available collection database [15] was beyond the scope of the present study but would be valuable.
Invasion success is often linked with a rapid rate of spread [58], and spread rate has been shown to correlate with factors such as the availability of vacant ecological niches [59].The shallow-water piscivore niche of B. belizanus has other members [28], although the fact that the species is immediately piscivorous after birth [20] and is resident year-round in such systems generally differentiates it from other piscivores.The niche (described as "piscivorous in open water habitat in the middle and high water column, principally using pectoral fins for motility") is recognized as being quite distinct from other fishes in the native range [60].The spread rate of B. belizanus estimated in the present study is intermediate to those of other invasive poeciliids.Invasive mosquitofishes (Gambusia affinis and G. holbrooki) spread much more quickly than B. belizanus, with rates ranging from ~38 km year −1 in the Colorado River watershed [59] to several hundred meters per day in small-scale studies [61,62].These rapid rates have been attributed to relatively high boldness of invasive mosquitofishes that contributes to their classification as some of the world's worst invasive species [58].In contrast, B. belizanus rate of spread is greater than the 1.6-km year −1 spread of Trinidadian Guppy Poecilia reticulata [59], another invasive poeciliid known for negative ecological effects [63], and similar to the 7.9 km year −1 spread of Sailfin Molly P. latipinna [59].Among non-poeciliids, the B. belizanus spread rate is comparable to that of Round Goby Neogobius melanostomus, for which estimates range from 0.5 to 17 km year −1 in the Great Lakes and Danube River [64][65][66].Some larger invasive fishes of particular concern exhibit very rapid spread rates, e.g., Nile Perch Lates niloticus (150 km 6 month −1 to 50 km week −1 ; [67], Bighead Carp (1.7-2.4 km day −1 ; [68]), and Northern Snakehead Channa argus (18 km month −1 ; [69]).Mayan Cichlid Cichlasoma urophthalmus spread from south to east-central Florida was estimated to be 15 km year −1 [70]; this species has been shown to negatively affect native fishes [71], with a range that includes Tampa Bay and is expanding [56,72].Depletion of prey in a given area by B. belizanus may increase the likelihood of range expansion [19], as may agonistic intraspecific interactions such as territoriality [73] and avoidance of cannibalism [14].
The disjunct distribution of the western population of B. belizanus in Tampa Bay relative to the eastern population is not likely to have resulted from swimming because over 15 km of deeper, open water separates the two populations and the species has not been collected in offshore habitat or in the many km of shoreline around the bay between the two populations (Figure 1).Data are lacking to determine when the species first occurred in western Tampa Bay, but potential mechanisms facilitating occurrence include aquarium release [74], bait bucket release [65,75], or bird carry-off [11,76].A single inseminated female B. belizanus released into a water body may be sufficient to establish a population because of large potential brood sizes (>300; [20]) and no evidence for deleterious phenotypic effects in poeciliid populations founded by a single individual [63].
Given the ability of B. belizanus to spread at up to 13 km year −1 , why then after over 20 years in Tampa Bay does the main, eastern population appear to be confined to a region that is a crow's-fly distance of ~30 km from north to south?Movement rates of fishes can be influenced by a number of factors, including the presence of suitable habitat [77].The B. belizanus range on the eastern fringe of Tampa Bay is bounded by large port facilities to the north and south (Figure 1) which are characterized by dredged ship channels and several km of extensive seawalls creating deep-water habitat at the shoreline.As shown in the habitat predictability analysis, the probability of presence of B. belizanus decreased with increasing depth.On the basis of the present study's results, it is hypothesized that the deep port seawall habitats with little cover could function as "soft" ecological barriers (sensu Bond and Lake [78]) to B. belizanus because of avoidance of such areas as a result of habitat preference [77] or because predation is high in the deeper water and no individuals survive [79].Between the two ports, the shoreline possesses a considerable amount of shallow-depth habitat including natural shorelines, mosquito ditches, and inland county drainage ditches, which probably facilitated spread of B. belizanus.The small western population may also be constrained by deep-water, seawalled habitats in the harbor adjacent to the creeks where the species was found in the present study.Spread of B. belizanus through the Everglades and south Florida was essentially unhindered by such potential soft barriers, at least to the south and west, and likely occurred through the extensive system of agricultural canals that include shallow nearshore areas.Thus, the species has reached the southwest coast to the first major urban area (near the city of Naples) but has not spread farther north (Figure 1) to other estuaries [80,81].The BRT habitat predictability analysis indicated that B. belizanus does occur along hardened shorelines (such as seawalls) in shallower water, albeit at somewhat lower probability of occurrence than other shoreline types (Figure 6), emphasizing the hypothesis that the deep-water aspect of seawalls in ports and harbors may be more important than the shoreline type.

Habitat Predictability
The present study displayed a common feature of species distribution modeling: the performance of the BRT models fitted to training data was higher than the cross-validation BRT models, which in turn was higher than the test data prediction model performance [7,46].Such patterns arise because prediction to test data is more difficult than simple interpolations to the sampled (training) regions, as a result of spatial habitat heterogeneity and varying ecological relationships [7].The BRT model of B. belizanus probability of presence based on the fixed subset of data had some value in its transferability, per the >0.6 AUC criterion of Huang and Frimpong [7], but the comparison to the BRT model based on the randomly selected subset of training data demonstrated that some of the apparent habitat relationships were idiosyncratic as opposed to general [8].The present study's emphasis on using prediction as a means of demonstrating ecological understanding [9] drove the need to consider model transferability not just in terms of cross-validation and prediction to test data, but also the consistency of observed species-habitat relationships through the comparison of the fixed and randomly selected subsets of data.Conducting 1000 iterations of the fixed and randomly selected BRT models allowed generation of 95% intervals, a novel approach that facilitated examination of the consistency of the observed relationships in order to evaluate which relationships appeared ecologically reasonable and therefore the predictability of habitat use.This approach is readily applicable to studies of other species where it is desirable to provide additional validation of observed species-habitat relationships.
Challenges in predictability of B. belizanus presence (Figure 4) may have reflected interannual differences in recruitment [37]-most waterbodies were sampled in 2006, but two were sampled in 2007-as well as differences between waterbodies in the equilibrium status of B. belizanus spread [82].In addition, B. belizanus is comparatively rare, even in the waterbodies with the highest percentage presence (Table 3).Rarity can hinder prediction of species-habitat relationships as a result of overfitting few presences with many covariates, suggesting that exploration of techniques such as ensemble-based modeling would be worthwhile in the future [83].Occupancy modeling based on repeated sampling at each site to account for imperfect detection may also provide additional insight [84].Initial colonization of new waterbodies can result in high abundance of B. belizanus, with numbers declining thereafter [85].Relatively low abundance generally appears to be typical of the species, even in the native range, and may reflect a conserved life history trait [37] which could be a function of some of the same factors influencing spread, e.g., territoriality [73], cannibalism [14], and resource depletion [19].
Dissolved oxygen was the second most important (after Julian day) of the habitat variables predicting B. belizanus presence in the studied estuarine waterbodies of Tampa Bay, although there was a considerable difference in the relationship between the BRT models based on the fixed and random subsets of data (Figure 6).The fixed data model showed occurrence was high at hypoxic conditions, reflecting a short-term event in Cockroach Creek during September 2007.Similar to other cyprinodontoids, B. belizanus may employ aquatic surface respiration in hypoxic conditions [86,87], possibly allowing the species to occupy such areas in relatively high numbers and prey successfully on species less adapted to low dissolved oxygen [88].Alternatively, aquatic surface respiration may have allowed B. belizanus to remain within low dissolved oxygen habitat, albeit with greater susceptibility to capture during sampling, akin to greater susceptibility to predation observed in other fishes exhibiting aquatic surface respiration [88].Although B. belizanus appears sufficiently physiologically flexible to cope with low dissolved oxygen, as observed in other invasive fishes in Florida [89], the relationships suggested in the present study were somewhat idiosyncratic rather than general and useful for prediction.
The greater probability of presence of B. belizanus with declining pH was consistent between fixed and random subset BRT models (Figure 6), and there are both similarities and differences with the patterns observed in previous studies.In freshwater habitat of south Florida's Everglades National Park, Kerfoot et al. [37] found a positive relationship between pH (range: ~7.4-8.0) and B. belizanus density at one site, a weak positive relationship at another site (pH range: ~7.2-7.4), and negative relationships at two other sites (pH ranges: ~7.3-7.6 and ~7.6-8.1).Multispecies ordination analyses within the native range do not demonstrate strong or consistent relationships between pH and B. belizanus populations [33][34][35][36].As noted by Kerfoot et al. [37], changes in pH are known to directly affect ion uptake across gills and internal acid-base balance in fishes, ultimately affecting survival rates.For example, low concentration of calcium ions was found to negatively affect the predation functional response of invasive N. melanogobius [90].Given potential future changes in pH as a result of climate change [91], future laboratory studies on B. belizanus could assess whether physiology and behavior (e.g., predation) are significantly affected by pH in a manner consistent with the lower probability of presence at higher pH observed in the present study.The Environmental Matching Hypothesis predicts that the impacts of an invader decline as habitat conditions move further from its physiological optimum [90].The present study found that B. belizanus is widely distributed across the salinity gradient, with evidence for greatest presence in polyhaline conditions (Figure 6).The observed effects of B. belizanus on small prey fish abundance in Tampa Bay estuarine waterbodies were consistent over a broad salinity gradient [24].The present study contrasts with that of Kerfoot et al. [37], who found decreasing B. belizanus density with increasing salinity in south Florida red mangrove habitats; a negative relationship with salinity has also been found in the native range [33].Water temperature in the coldest month is 2-3 • C cooler in the Tampa Bay watershed than south Florida [92], so that occupation of higher salinity by B. belizanus could ameliorate low-temperature effects.Such an effect was found for Tilapia mossambica (now Oreochromis mossambicus, a euryhaline freshwater species), which is restricted to estuaries near the southern limit of its range in Africa, possibly in order to maintain near-normal sodium and chloride ion concentrations during low winter temperatures [93].A similar mechanism was not found for C. urophthalmus, however [94].Nekton abundance often increases with increasing salinity because of the increasing contribution of marine-derived juveniles [43,95] and this may provide a greater potential pool of prey resources for B. belizanus at higher salinity.The estuarine focus of the sampling in the present study precludes any firm conclusions about the distribution or abundance of B. belizanus in the permanent, non-tidal freshwater reaches upstream of the study area; expansion of sampling in this habitat using similar seining methods would clarify the significance of the present study's findings with respect to salinity.
Lower lethal temperature (~9.5 • C) for B. belizanus is surprisingly low and rarely experienced in Tampa Bay [96]; feeding is, however, reduced (or ceases) at ~11-16 • C [92,97], which are more commonly observed winter temperatures in the area.During an 80-year cold front in south Florida, water temperature in some areas dropped below the laboratory-determined B. belizanus mean lethal temperature, yet the species persisted [98].It would be informative to repeat Shafland and Pestrak's [92] and Kerfoot's [97] experiments, which were undertaken in freshwater aquaria, at a variety of salinity levels to assess potential ameliorating effects on low temperature effects.The present study found a consistent drop-off in B. belizanus probability of presence below about 19.5 • C (Figure 6), whereas Kerfoot et al. [37] found that density of B. belizanus in south Florida red mangroves increased as temperature decreased, although the lowest temperature in their dataset was 23 • C. Long-term temperature increases may favor the species' increased expansion [38], as noted for C. uropthalmus [71], producing effects on native species that may not yet have been manifested [99].

Management Implications
The range of shoreline habitats and physicochemical conditions at which B. belizanus was present in this study suggests that within certain parameters (i.e., nearshore, shallow-water habitat), the species has generally broad environmental tolerance, a key trait of successful invasive fishes [6].The consistent importance of Julian day in predicting probability of presence suggests that future studies aiming to assess the occurrence of B. belizanus in Tampa Bay estuarine waterbodies should be focused on the period from September to December/January in order to maximize the likelihood of capture (Figure 6).The present study suggests that B. belizanus in estuarine areas of Tampa Bay spread at a modest rate compared to the most notable invasive fishes and that its distribution may be constrained by soft ecological barriers, i.e., deep seawall habitat.The recent classification as moderately invasive [23] appears appropriate.Nevertheless, strong evidence for negative effects on small-bodied fishes [24] and potential competition with C. undecimalis [15] indicates that consideration of possible management actions is warranted.Under current climate conditions, B. belizanus would be predicted to be able to occur 150-200 km farther north than Tampa Bay based on lower lethal temperature [92,97].However, the present study's hypothesis that spread of B. belizanus is restricted by deep, anthropogenically altered habitats leads to the conclusion that movement north may be very much dependent on human-facilitated pathways, recently referred to as "jump dispersal" [100], for which the main vectors are releases from aquaria and bait buckets [75].Internet-based sales of ornamental fish allow species to be spread easily within the U.S., suggesting that education of hobbyists and retailers is paramount [74] and particularly important in the case of B. belizanus, given that aggression by piscivores is one reason for release of fish from aquaria [101].
Regarding bait bucket releases, B. belizanus is not raised commercially as bait, but releases could occur as surplus fish from anglers cast-netting [102] and catching the species along with other small-bodied fishes.This is particularly true within the eastern portion of the B. belizanus Tampa Bay range, which is notably good habitat for sub-adults and adults of the primary inshore sportfish C. undecimalis [103].One possible means of limiting spread by this vector would be to produce simple educational materials for anglers (e.g., posters at boat ramps) illustrating the current range of B. belizanus, emphasizing the need to not release it beyond this range, and giving warnings about its potential effects through competition with prized sportfishes that anglers may be targeting.Such materials could be combined with similar warnings for other spreading invasive fishes, such as C. uropthalmus [56].Emphasis on not releasing B. belizanus from aquaria or bait buckets is also relevant for freshwater habitat because of the similarity in prey capture and potential competition with juveniles of another prized sportfish, Largemouth Bass Micropterus salmoides [104].
Restoration is being undertaken in the Tampa Bay watershed with the goal of approximately restoring the balance of estuarine habitat observed ca.1950 [41].Projections of habitat change as a result of sea level rise suggest that mangrove habitat will increase its percentage of coastal habitat coverage from ~74% to 85-89% by 2100, with a concomitant decrease in salt marsh from 25% to 10-14% [41].The habitat predictability analysis demonstrated that B. belizanus is sufficiently flexible in terms of shoreline habitat occupation (Figure 6) that such changes would not be expected to result in considerable effects to the species.With a spread rate of several km year −1 , it would be expected that B. belizanus would colonize restored areas relatively quickly, particularly because some of the major restored sites are directly upstream of two of the waterbodies in the present study, Andrews Creek and Cockroach Creek [105].This could disrupt typical successional fish assemblage dynamics that are observed in restored areas [106][107][108], as evidenced by effects on small-bodied resident fishes elsewhere in Tampa Bay [24] and potential competition for prey with transients using restored areas as a nursery, in particular C. undecimalis [15].

FishesFigure 1 .
Figure 1.Presence of Belonesox belizanus in estuarine waterbodies of Tampa Bay, 1989-2007/2008, amended with permission from [24], with inset showing collections from elsewhere in Florida [15].For main map, red dots represent presence (with a red star indicating the first occurrence, 1994), green dots represent absence.QC = Question Mark Creek; RW = Riverview West Creek; DC = Dogleg Creek; MBC = Mill Bayou Canal; SCS = Sun City Slough.Total n = 16,113 shoreline seines (see Supplementary Material).Underlined waterbodies were included in the boosted regression tree modeling of habitat predictability.Collections from elsewhere in Florida represent clustered records ranging from one specimen (white circle) to 20 or more specimens (dark blue circle).

Figure 1 .Fishes
Figure 1.Presence of Belonesox belizanus in estuarine waterbodies of Tampa Bay, 1989-2007/2008, amended with permission from [24], with inset showing collections from elsewhere in Florida [15].For main map, red dots represent presence (with a red star indicating the first occurrence, 1994), green dots represent absence.QC = Question Mark Creek; RW = Riverview West Creek; DC = Dogleg Creek; MBC = Mill Bayou Canal; SCS = Sun City Slough.Total n = 16,113 shoreline seines (see Supplementary Material).Underlined waterbodies were included in the boosted regression tree modeling of habitat predictability.Collections from elsewhere in Florida represent clustered records ranging from one specimen (white circle) to 20 or more specimens (dark blue circle).

Figure 2 .
Figure 2. Presence (red dots with last two digits of year annotation) of B. belizanus in the lower Little Manatee River sampling area (km 0-13) from regular shoreline seine monitoring, 1997-2007.The first six presences are highlighted, along with the most upstream presence.The broken yellow line divides the upstream and downstream zones of monthly sampling.The inset graph shows rate of upstream spread from the first to sixth occurrence in 1997/1998.

Figure 2 .
Figure 2. Presence (red dots with last two digits of year annotation) of B. belizanus in the lower Little Manatee River sampling area (km 0-13) from regular shoreline seine monitoring, 1997-2007.The first six presences are highlighted, along with the most upstream presence.The broken yellow line divides the upstream and downstream zones of monthly sampling.The inset graph shows rate of upstream spread from the first to sixth occurrence in 1997/1998.

Fishes 2017, 2 , 23 Figure 3 .
Figure 3. Predictive performance of boosted regression tree models for B. belizanus presence in Tampa Bay estuarine waterbodies for randomly selected (|---|) and fixed (|---|) subsets of 935 9.1-m shoreline seine samples, expressed as 95% intervals of the area under the receiver operating characteristic curve from 1000 model iterations.Randomly selected and fixed subset models included training data (n = 695 samples) and prediction data (n = 240 samples).

Figure 3 .
Figure 3. Predictive performance of boosted regression tree models for B. belizanus presence in Tampa Bay estuarine waterbodies for randomly selected (|---|) and fixed (|--|) subsets of 935 9.1-m shoreline seine samples, expressed as 95% intervals of the area under the receiver operating characteristic curve from 1000 model iterations.Randomly selected and fixed subset models included training data (n = 695 samples) and prediction data (n = 240 samples).

Fishes 2017, 2 , 23 Figure 3 .
Figure 3. Predictive performance of boosted regression tree models for B. belizanus presence in Tampa Bay estuarine waterbodies for randomly selected (|---|) and fixed (|---|) subsets of 935 9.1-m shoreline seine samples, expressed as 95% intervals of the area under the receiver operating characteristic curve from 1000 model iterations.Randomly selected and fixed subset models included training data (n = 695 samples) and prediction data (n = 240 samples).

Table 1 .
Tampa Bay estuarine waterbodies sampled to assess Belonesox belizanus habitat predictability, including mean and range of continuous habitat predictors.

Table 2 .
Shoreline and substrate type in Tampa Bay estuarine waterbodies sampled to assess B. belizanus habitat predictability.

Table 3 .
Abundance and percentage presence of B. belizanus in Tampa Bay estuarine waterbodies sampled to assess habitat predictability.