Natural abiotic factors more than anthropogenic perturbation shape the invasion of Eastern Mosquitofish (Gambusia holbrooki)

Eastern Mosquitofish (Gambusia holbrooki) is an invasive and globally widespread species that is considered highly tolerant. We used species distribution models (SDMs) to assess factors, including the role of anthropogenic perturbation, that mediate its invasion on a regional scale. A better understanding of the important large-scale factors may help us identify future areas of concern and potential avenues for control. We built SDMs from presence records and randomly selected pseudo-absences of mosquitofish on a 10- × 10-km grid. We used 10 modeling techniques implemented in the biomod2 software and ensemble forecasts. Final models contained 12 environmental predictors, including natural environmental factors (elevation, slope, topographic index, precipitation, accumulated flow, temperature mean and range) and anthropogenic perturbation indicators (population density, urban and agricultural land uses, number of local and upstream dams). Elevation, temperature, and accumulated flow most strongly influenced mosquitofish distribution, and mosquitofish were found more frequently in downstream, warmer waters. Anthropogenic features, except the number of upstream dams, were poorer predictors of mosquitofish presence than were natural environmental factors. The best models suggested that mosquitofish are more likely to occur in areas with more dams upstream, but removing this predictor did not strongly affect model results. Restoration efforts or modifications to anthropogenic features appear unlikely to alter mosquitofish distribution patterns, highlighting the importance of preventing introductions to new areas. Mosquitofish have been extensively documented in the Iberian Peninsula, but consensus methods suggest many additional suitable areas from which records were not found. Thus, this highly invasive species is or may become much more widely distributed than current observations in the region.

Eastern and Western Mosquitofish (Gambusia holbrooki and the closely related Gambusia affinis) are among the world's 100 most invasive species (Lowe et al. 2000). They have expanded from their original range in part of the USA and Mexico to become among the most prolific, widespread freshwater fishes in the world (García-Berthou et al. 2005, Pyke 2008). Mosquitofish can survive a wide range of physical conditions and may reproduce very rapidly (Pyke 2008). Their invasions are marked by ecological changes. These 2 poeciliids have been implicated globally in severe amphibian declines (Gamradt andKats 1996, Baber andBabbitt 2004), local extirpations of native fishes (Alcaraz et al. 2008), and trophic cascade effects (Hurlbert et al. 1972, Hurlbert andMulla 1981). Predation on zooplankton and insects by mosquitofish has caused cascading changes to phytoplankton and primary production with effects on light, nutrients, temperature, pH, and dissolved O 2 (Hurlbert et al. 1972, Hurlbert andMulla 1981). Mosquitofish have been intentionally introduced across the globe, in part because of the frequent misconception that they have a restricted diet and are more efficient predators of mosquitoes than native species (see Pyke 2008 and references therein). Once introduced, we assume that abiotic factors, such as temperature and flow, shaped their ability to become established and to invade a landscape by changing environmental resistance, as predicted by Moyle and Light (1996) based on observations of invasive fishes in California.
Invasibility, the susceptibility of environments to invasion by introduced species, such as mosquitofish, is im-portant in understanding current and future species distributions and in explaining the factors that enhance and limit invasions (Lonsdale 1999, Alpert et al. 2000. Temperature, elevation, and water flow strongly affect mosquitofish survival and reproduction at local scales, so we might expect these factors to influence the range and extent of mosquitofish on a landscape scale. For example, temperature is a determinant of the latitudinal extent of mosquitofish through its influence on reproductive rates, overwintering success, and biotic interactions, such as competition with the native Iberian Toothcarp Aphanius iberus (Benejam et al. 2009. Higher elevation typically is associated with lower water temperatures, steeper stream habitats, and reduced occurrences of the backwater and side-channel habitats frequently occupied by mosquitofish. The influence of flow on its preferred habitats (Meffe 1984) explains why mosquitofish often are found in areas with low or regulated flow conditions, where warm backwaters can support planktivorous feeding and productivity. We hypothesized that elevation, temperature, and water flow are critical factors in mosquitofish distribution across a landscape.
Human activities, such as land use and construction of dams and reservoirs, and their ecological consequences, can dramatically alter invasibility, but they often are neglected in considerations of invasions over broad extents. This neglect seems counterintuitive because human effects themselves may be easier to quantify over broad extents than their specific ecological consequences. On fine scales, fish assemblages commonly are used in biological indices as metrics of stream health and anthropogenic impacts on stream quality (Hughes et al. 1998), which can influence species distributions (e.g., Vila-Gispert et al. 2002). We expect these relationships to be direct through changes in habitat suitability and indirect through alteration of the biotic community structure and effects on other species.
Mosquitofish can thrive in streams with altered flow regimes and are rare in cold headwater streams , in part because of temperature, but also because they are limnophilic species that do not withstand strong water flows well (Meffe 1984). Mosquitofish are highly tolerant and can thrive in agricultural ditches, sloughs, and urban streams. Thus, these areas may provide additional habitat opportunities (Cabral and Marques 1999, Baber et al. 2002, Northington and Hershey 2006, Simon and Travis 2011. In Mediterranean systems, large dams regulate flow regimes and reduce extreme flow events (Bunn and Arthington 2002). Dams limit high flows and create reservoirs that can provide a constant source of downstream propagules if they are inhabited by mosquitofish. Freshwater exotic fishes often are introduced through reservoirs, from which they spread to rivers (Rahel 2002). We hypothesized that hydrological alterations, such dams and reservoirs, would favor mosquitofish invasion and that mosquitofish were more likely to be introduced where human activities were more apparent, such as in areas with urban or agricultural development.
We used species distribution models (SDMs) to assess how natural environmental factors and anthropogenic disturbances mediate invasion by mosquitofish, and we used the Iberian Peninsula as a case study. Our objectives were to: 1) understand the factors that mediate mosquitofish invasion, and 2) test the role of anthropogenic disturbances, such as hydrological alteration, in the invasion process. We hypothesized positive relationships between mosquitofish presence and anthropogenically driven hydrological alteration and agricultural land use. As far as we know, we are the first to present an SDM for a mosquitofish species. Mosquitofish have been generally neglected in similar models, except by Lauzeral et al. (2011) who noted some climate niche shifts. Our SDM is one of the first in which anthropogenic landscape features were considered for a freshwater species (see also Markovic et al. 2012, Liang et al. 2013, Maloney et al. 2013. We hope to improve our understanding of factors mediating mosquitofish invasion and to identify areas where mosquitofish may already be present or is likely to invade.

Data compilation
We obtained 919 unique records from >1000 documented occurrences of mosquitofish in the Iberian Peninsula, an area of ∼582,000 km 2 including Spain and Portugal, by reviewing research articles, atlases, and databases (published and unpublished), to build a map with documented presence data (e.g., Doadrio 2002, Ribeiro et al. 2007; see Appendix S1 for additional sources). We considered all mosquitofish (Gambusia spp.) records to be the Eastern Mosquitofish (G. holbrooki) because recent genetic studies have identified only G. holbrooki in the Iberian Peninsula and only G. holbrooki is known to have been introduced (Vidal et al. 2010). This distinction is important because grouping of G. affinis and G. holbrooki as subspecies of G. affinis until around 1988 resulted in uncertainty in their respective distributions and inconsistent naming (Wooten et al. 1988, Rauchenberger 1989, Robins et al. 1991, and some references still cite G. affinis. We used a resolution of 10 km UTM (Universal Transverse Mercator; i.e., 100 km 2 ) because it was the most widely used resolution in the records and could be obtained for all data. Our final grid covered the entire Iberian Peninsula and contained 6138 unique 10 km UTM grid cells. Mosquitofish were recorded in ∼15% of the cells (Fig. 1A).
We included variables that have been frequently used in SDMs of freshwater fish, such as elevation, topographic index, flow accumulation, slope, land cover, precipitation, annual thermal range (the difference between annual minimum and maximum temperatures), and temperature, as predictors of mosquitofish presence (e.g., McNyset 2005, Chen et al. 2007, DeVaney et al. 2009, Grenouillet et al. 2011. We compiled environmental data layers from online databases in ArcGIS (version 10; Environmental Sys-tems Research Institute, Redlands, California) and did subsequent calculations with its Spatial Analyst toolbox. We considered natural descriptors thought to be related to the establishment, probability, and abundance of mosquitofish (e.g., temperature [ Fig ber of dams [ Fig. 1D], landuse alterations) thought to affect the likelihood of introduction, establishment, or abundance and, thus, likelihood of occupancy. Factors affecting abundance would be expected to influence survival and reproduction and might suggest factors that, at their extremes, are highly suitable or unsuitable for this species.
When environmental data were originally of finer grain than our species data, we used the mean of values within each 10-× 10-km grid cell. As recommended by Dormann et al. (2013), we removed highly correlated variables (| r | ≥ 0.7), keeping the variable that we expected to contain the most relevant information (e.g., mean annual temperature was retained rather than minimum temperature of the coldest month or maximum temperature of the warmest month when all 3 were correlated). Removing highly correlated variables reduces the resulting variance of regression parameters, improves processing time, reduces errors during processing, and prevents possible misinterpretation of results Thuiller 2005, Elith andLeathwick 2009a). However, if predictor variables have independent ecological relevance, they should be retained. According to a recent, comprehensive review (Dormann et al. 2013), regression-type approaches (e.g., generalized linear models) and machine-learning methods work reasonably well under moderate collinearity and are not outperformed by techniques specifically designed to deal with collinearity (e.g., partial least squares). Therefore, neither elevation nor mean annual temperature were eliminated for the SDMs because, although correlated, both have independent effects on mosquitofish populations (Carmona-Catot et al. 2011). In the Iberian Peninsula, this correlation is presumably exacerbated by topography, including the location of the Pyrenees and Cantabric mountain ranges to the north. Our original environmental matrix had 23 variables, but only 12 were used (Table 1) after elimination of highly correlated predictors.

Modeling
For fitting species distribution models, we used the BIOMOD computation framework (Thuiller 2003) as implemented in the biomod2 package (Thuiller et al. 2013) in R (version 3.0.3; R Project for Statistical Computing, Vienna, Austria). We used generalized linear models (GLM), generalized boosting models (also known as boosted regression trees, BRT), classification tree analysis (CTA), mixture discriminant analysis (MDA), random forests (RF), artificial neural networks (ANN), generalized additive models (GAM), surface range envelopes (SRE), multiple adaptive regression splines (MARS), and Maxent (see Olden et al. 2008, Elith and Leathwick 2009b, Franklin and Miller 2009, Dormann et al. 2013 for reviews of these methods). On average, the different modeling methods yield good agreement between observed and predicted distributions, but relative performance of different techniques is idiosyncratic across species (Thuiller 2003). The BIOMOD approach implements several widely used regression techniques (e.g., GLM, GAM, and MARS). It also includes a number of more modern machine learning methods, such as RF and BRT, which offer greater predictive performance for nonlinear relationships, compared to the more traditional methods (Olden et al. 2008, Elith andLeathwick 2009b). Combining the results of the different techniques with ensemble forecasting often yields more robust predictions New 2007, Thuiller et al. 2009).
The mosquitofish presence data mostly originated from planned electrofishing quantitative surveys with accurate location records, but we treated the data as presence-only and built pseudo-absences because fish species absence is very difficult to ascertain accurately and probably is subject to a number of potential biases (e.g., contrasting surveying effort; effects of stream size, water conductivity, and flow on detectability; Barbet-Massin et al. 2012). For tests run in biomod2, we used 1000 randomly selected pseudoabsences, with no repetition of selection, for 2-run evaluations. These values were estimated based on the areal coverage recommendations provided in Barbet-Massin et al. (2012). We split our data into 70/30 training and test data and randomly selected pseudo-absences. We compared the results to those obtained with repetition and ½ the number of randomly selected pseudo-absences and noted no obvious deviations. As such, only results with the initial settings are presented below.
We evaluated model accuracy (predictive power) of the different techniques using the following criteria: the area under the receiver operating characteristic (ROC) curve (AUC), Cohen's kappa (κ), the true skill statistic (TSS), and sensitivity (the proportion of correctly predicted presences) and specificity (the proportion of correctly predicted absences) for each (Allouche et al. 2006, Thuiller et al. 2009), where the binary threshold was defined as the one that maximized the evaluation metric value. As recommended by Allouche et al. (2006), we focused on TSS over κ because the former does not depend on prevalence (the proportion of sites in which the species was recorded as present). We compared the variable contributions from all models to determine the most influential environmental factors. We noted the relative magnitude and order of variable contributions, and we paid special attention to the models with best accuracy. Last, we computed an ensemble forecast built using models with a TSS score greater than 0.7 (RF and Maxent) weighted by their evaluation scores. We also built ensemble forecasts with the same threshold for AUC and κ, only to evaluate these metrics on ensemble performance.

RESULTS
The model outputs confirmed our hypotheses based on the previous knowledge of mosquitofish ecology and predicted an increasing probability of mosquitofish presence at lower altitudes with lower gradients and higher mean temperatures. The model also predicted a trend toward increasing probability of mosquitofish presence with reduced precipitation and increased thermal range. With the exception of cumulative upstream dams, anthropogenic factors had little to no effect on the probability of mosquitofish presence.
The SDM results, including individual and ensemble models and especially the best-performing model (RF), suggested that mosquitofish could occupy a wider extent than has been documented in the Iberian Peninsula (Fig. 1E, F). The RF model predicted that >½ of the Iberian Peninsula, including much of the southwest from where few records exist, appears highly suitable for mosquitofish (Fig. 1E).
The importance of each variable in all 10 model calibrations was more consistent for variables that were either very strong or very weak predictors (Table 2). Elevation was the most influential variable in all models explaining mosquitofish distribution, but its contribution to each model varied. Intermediate variables were not ranked consistently across tests. The 2 nd -most influential variable varied by test and included mean annual temperature, thermal range, slope, annual precipitation, and flow accumulation. All of those variables, except flow accumulation, appeared in the top 5 variables by contribution for ≥½ of the tests. Among environmental drivers, topographic index was the poorest predictor and was not notably important in its contribution to any model. At best, it contributed just >3%. Flow accumulation also generally contributed little to model building, except in the case of ANN. The other environmental predictors (mean annual temperature, annual precipitation, thermal range, slope, and elevation) were all stronger contributors to ≥½ of the calibrated models (Table 2).
Anthropogenic variables generally were poor predictors compared to abiotic drivers, such as elevation, temperature, and thermal range ( Table 2). The cumulative-upstreamdams metric was the only anthropogenic variable ranked in the top 3 predictors (RF, Maxent, and CTA models). Models run without cumulative upstream dams revealed that it had little effect on model performance. Of the landuse measures, agriculture appeared to be a better predictor than urban use. However, % agricultural land use was in the top 5 variables for <½ of the tests and always had a low contribution (<10%). Number of local dams and % local urban land use were always poor predictors of mosquitofish presence and were not included in the top 5 variable contributions for any of the 10 tests. Human population density and flow accumulation also were poor predictors, contributing only to ANN and SRE results, which appeared to deviate from the general patterns of variable rankings observed across other tests.
RF and Maxent, followed by BRT, were the best models (based on TSS and AUC criteria; Table 3). RF and Maxent were the only models with TSS scores >0.7. In addition to their high total scores, RF and Maxent were consistently among the most sensitive and the most specific models judged by all evaluation metrics. In contrast, the ensemble model performed poorly relative to best performing models. Four of the 5 most important predictors were shared by RF, Maxent, and BRT: elevation, slope, total number of upstream dams, and mean annual temperature (Table 2). Thermal range was important in the RF model, whereas annual precipitation was important in the Maxent model. In the RF model, elevation was the strongest contributor (20.7%), followed closely by mean annual temperature (17.5%). Slope (8.3%), total number of upstream dams (7.9%), thermal range (7.3%), flow accumulation (6.2%), and annual precipitation (6%) all contributed similarly, but human population density, land use (urban and agriculture), local dams, and topographic index contributed little.

DISCUSSION
The Iberian Peninsula presents a unique opportunity to study mosquitofish invasions because it was the site of the first mosquitofish introduction in Europe (in 1921). Mosquitofish quickly invaded the Mediterranean basin after this intentional introduction. By 1948, they appear to Peninsula also has a landscape over which diverse anthropogenic effects can be assessed. In the region, extensive modification and flow regulation have resulted in documented changes to habitat and water-quality characteristics, producing aquatic systems that in some cases are dramatically different from pristine conditions. Invasibility varies with both natural and anthropogenic disturbances. Freshwater ecosystems are regarded as inherently more susceptible than other ecosystems to invasion (Shea and Chesson 2002), and increases in invasibility probably have resulted from anthropogenic alterations, such as those to water quality, flow regimes, and substrate composition (Mooney and Hobbs 2000, Marchetti and Moyle 2001, Marchetti et al. 2004, Lapointe et al. 2012). The consequences of increased invasibility might be especially important in Mediterranean climates, where extensive human impacts and exotic fish introductions have occurred in regions originally supporting a high proportion of endemic freshwater species (Marr et al. 2010, Ribeiro andLeunda 2012). The Eastern Mosquitofish has been implicated regionally in declines of the endemic, endangered cyprinodontiforms Valencia hispanica and Aphanius spp. and top-down effects in shallow Mediterranean lakes (Blanco et al. 2004, Caiola and de Sostoa 2005, Alcaraz et al. 2008, reinforcing the importance of understanding its ecology and distribution in conservation and restoration efforts. Landscape-scale factors that appear to be important to mosquitofish distribution were consistent with our expectations. Mosquitofish prefer low-flow habitats (Meffe 1984); temperature plays a critical role in their survival and reproduction (Benejam et al. 2009) at local scales, and they are excluded from high latitudes (Pyke 2008). Therefore, our finding that elevation was the strongest and most consistent contributor to SDMs seems a natural outcome. We would expect elevation to represent gradients of available side-channel and backwater habitats, temperature, and water velocity and that all of these variables (fewer backwaters, lower temperatures, and water velocity) would make life at high elevations less than ideal for this fish. The consistency across tests with which environmental drivers were more influential than anthropogenic factors leaves us confident of the importance of these broad environmental factors to mosquitofish.
Our results agree with previous reports that machine learning methods such as RF, Maxent, and BRT are among the best-performing (e.g., Elith et al. 2006, Franklin andMiller 2009). The models with these techniques had very good accuracy, sensitivity, and specificity as calculated based on their AUC, κ, and TSS scores. These models strongly support environmental gradients as critical to mosquitofish invasion success. The importance of environmental gradients also is broadly consistent with recent work by Strecker and Olden (2014), who found that phylogenetic clustering in nonnative species was related to natural environmental conditions.
For mosquitofish, human impacts appear to be less important than broad-scale environmental gradients. For example, agricultural land use was not an important pre- Table 3. Performance of the different modeling techniques for the species distribution models of Eastern Mosquitofish in the Iberian Peninsula. GLM = generalized linear models, BRT = boosted regression trees, CTA = classification tree analysis, MDA = mixture discriminant analysis, RF = random forests, ANN = artificial neural networks, GAM = generalized additive models, SRE = surface range envelopes, MARS = multiple adaptive regression splines, AUC = area under receiver operating characteristic curve, κ = Cohen's kappa, TSS = true skill statistic, sensitivity = the proportion of correctly predicted presences, specificity = proportion of correctly predicted absences. dictor, possibly because the models produce probabilistic maps of spatial occupation without predicting or reflecting mosquitofish densities in occupied areas. Nevertheless, this finding suggests that changing anthropogenic landscape features may yield only limited gains in slowing or preventing the spread of mosquitofish, at least in the Iberian Peninsula. The number of upstream dams is an indicator of flow regulation. In many systems, especially Mediterranean waterways prone to flash flooding, dams dampen the peaks and troughs characteristic of natural flow regimes (Lytle andPoff 2004, Poff et al. 2007). Based on the habitat preferences of mosquitofish, we expected the number of upstream dams to have a stronger effect on models of mosquitofish distribution. It is unlikely that mosquitofish would be able to maintain positions in systems where it would be regularly flushed out of headwaters. In addition, its competitive dominance over local species appears to be lower where salinity fluctuations are more extreme, i.e., in unregulated downstream reaches near estuarine and oceanic conditions (Alcaraz et al. 2008). One could reason that restoring the natural flow regime or removing dams might decrease the abundance and impact of mosquitofish. However, we expect that even if mosquitofish could be temporarily removed from an area, such action would be unlikely to eradicate it, in particular because it could recolonize from refuge habitats, such as ponds and floodplains.
We had reasoned that after its long history of introductions, i.e., a high propagule pressure, mosquitofish may have reached the apex of its invasion in the Iberian Peninsula, and we assumed it was maintaining an extent more or less at equilibrium with abiotic factors. Whether this equilibrium actually exists is unclear because the bestperforming models suggest that mosquitofish distribution is grossly underestimated and that its expected distribution is >4× what has been recorded. In general, the model results appear to identify tributaries to major river systems and the southwestern corner of the Iberian Peninsula as areas with high habitat suitability, suggesting that these areas might harbor undocumented populations or might already have been invaded by mosquitofish. Its long history of invasion would leave us to suspect that it is already present in at least some of these locations, and we would strongly recommend that researchers and management agencies check these areas for established populations. Mosquitofish are small and may be more easily overlooked than other invasive species. An alternative possibility is that our models are missing a landscape-scale factor that limits mosquitofish distributions in the areas where our current models predict them, leading to an overestimation of suitable area. Based on our ecological knowledge of the species, this possibility seems less likely than a lack of observation or previous introduction, especially because scattered records exist near these areas (Fig. 1D).
Our results suggest that natural abiotic factors are much more important than hydrological disturbance in regulating mosquitofish presence. Models of climate change predict that the region probably will experience warming, concentrated during the summer months, and markedly decreased precipitation with increased variability (Álvarez Cobelas et al. 2005). Estimates specific to inland waters suggest these changes will have stronger effects on endorheic and groundwater systems, alpine waters, and coastal wetlands than on other water bodies (de Castro et al. 2005, Iglesias et al. 2005. Less predictable flow events have the potential to flush mosquitofish downstream, but we think it unlikely that they would eliminate populations because hydrological regulation (in the form of dams and reservoirs) was less important than expected in our models. However, the potential for warming in uninvaded high-altitude areas is of great concern. Temperature and elevation appear to be extremely important limiting factors for the current distribution of mosquitofish in the region. Warmer temperatures and reduced flows could mean that new areas will become suitable for mosquitofish in the future, thereby aiding its invasion. The relative unimportance of anthropogenic hydrological alterations also suggests that restoration efforts might reduce mosquitofish abundance, but are unlikely to extirpate the species from an area. Thus, prevention of further introductions must become a priority. We urge extreme caution, especially in the areas where mosquitofish presence has been predicted but not observed. Given the importance of abiotic rather than anthropogenic factors as limitations on its distribution, once mosquitofish are introduced in a suitable region, there will be few realistic opportunities to make large changes that can limit its distribution.