Off to new shores: Climate niche expansion in invasive mosquitofish (Gambusia spp.)

Abstract Aim Formerly introduced for their presumed value in controlling mosquito‐borne diseases, the two mosquitofish Gambusia affinis and G. holbrooki (Poeciliidae) are now among the world's most widespread invasive alien species, negatively impacting aquatic ecosystems around the world. These inconspicuous freshwater fish are, once their presence is noticed, difficult to eradicate. It is, therefore, of utmost importance to assess their geographic potential and to identify their likely ability to persist under novel climatic conditions. Location Global. Methods We build species distribution models using occurrence data from the native and introduced distribution ranges to identify putative niche shifts and further ascertain the areas climatically suitable for the establishment and possible spread of mosquitofish. Results We found significant niche expansions into climatic regions outside their natural climatic conditions, emphasizing the importance of integrating climatic niches of both native and invasive ranges into projections. In particular, there was a marked shift toward tropical regions in Asia and a clear niche shift of European G. holbrooki. This ecological flexibility partly explains the massive success of the two species, and substantially increases the risk for further range expansion. We also showed that the potential for additional expansion resulting from climate change is enormous—especially in Europe. Main conclusions Despite the successful invasion history and ongoing range expansions, many countries still lack proper preventive measures. Thus, we urge policy makers to carefully evaluate the risk both mosquitofish pose to a particular area and to initiate appropriate management strategies.


| INTRODUC TI ON
Globalization with its massive global trade and long-distance transportation is leading to a steady increase in the number of biological invasions, affecting all taxonomic groups and all continents, with no sign of saturation (Seebens et al., 2017). The rapidly changing climate is further facilitating the spread and establishment of invasive alien species (IAS; Hulme, 2017). IAS represent a major threat to biodiversity, challenging conservation efforts and management of biological resources (Simberloff et al., 2013). Accordingly, the Convention on Biological Diversity's (CBD) Strategic Plan for Biodiversity demanded a substantial increase in efforts made to reduce the impact and spread of invasive species, with prioritizing global actions of management and control (Essl et al., 2020;McGeoch et al., 2016).
However, an essential prerequisite for the management and control of IAS is to understand the factors that determine the geographic distribution of a species and prevent it from spreading to other ecosystems.
The eastern mosquitofish, Gambusia holbrooki, and the closely related western mosquitofish, G. affinis, are one of the most successful freshwater IAS. They are native to the eastern and central United States, respectively, but have been introduced to every continent except Antarctica by aggressive introduction programs and their presumed value as mosquito control agents (Fryxell et al., in press;Krumholz, 1948;Pyke, 2008;Stockwell & Henkanaththegedara, 2011). Both species are tolerant toward anthropogenic disturbances (e.g., pesticides) and are capable of surviving a broad range of environmental conditions, as exemplified by tolerating salinities up to 41 ppt (Hubbs, 2000), temperatures between 0 and 40°C (Cherry et al., 1976;Lau et al., 2019), or oxygen contents ranging well into the hypoxic range (Cherry et al., 1976;Odum & Caldwell, 1955;Santi et al., 2020). These characteristics along with bearing live young contribute to their success as invasive species (Pyke, 2008;Walton et al., 2012). Collectively, G. holbrooki and G. affinis are among the most invasive fish worldwide and are currently considered as one of the 100 most detrimental IAS (Lowe et al., 2000). Their negative impact on local faunas stems partially from their often carnivorous feeding behavior (Pirroni et al., 2021;Pyke, 2008), and indigenous fish and amphibian larvae often rapidly decline after the introduction of mosquitofish (Barrier & Hicks, 1994;Morgan & Buttemer, 1996;Remon et al., 2016). Their dramatic effect for the local (often endemic) fauna has now been widely documented, especially in Australia (Arthington, 1989;Ivantsoff, 1999) and Europe (Alcaraz et al., 2008;Alcaraz & García-Berthou, 2007;Carmona-Catot et al., 2013;Rincon et al., 2002).
Despite this, Gambusia spp. together with other Poeciliid fishes (e.g., Poecilia reticulata) are still used as mosquito control agents in some parts of the world (Jayapriya & Shoba, 2014;Saleeza et al., 2014;Verma et al., 2016). From a conservation perspective it is, therefore, essential to identify regions where the-deliberate or accidental-introduction results in a high probability of establishment due to suitable (current and future) climatic conditions. Moreover, changing boundaries of already established populations also need to be robustly assessed. This is especially important, because prevention of spread (e.g., via control and public education) is more effective than trying to eradicate established populations (Fournier et al., 2019). A great tool for this type of assessment are ecological niche models (ENM; Guisan & Thuiller, 2005;Guisan & Zimmermann, 2000). ENMs are correlational techniques aimed to identify climatic regions where the species might find suitable conditions, based on current occurrence data and climatic information, and are broadly applied in the fields of biology, nature conservation, and biogeography (Elith et al., 2011). However, these models typically rely on the assumption that species retain their climatic niche in the exotic ranges (Peterson, 2011;Wiens et al., 2010). While this is true for many species, the full climatic potential is often not even fully realized due to sporadic introductions and a limited dispersal capacity (Pearson, 2007;Sillero, 2011). This is particularly pertinent for aquatic ecosystems where species mostly disperse within the river networks (Tonkin et al., 2018). Yet, some invasive species may even occur under climatic conditions that are outside the range of climatic values they inhabit within their native geographic distribution (Broennimann et al., 2007;Medley, 2010;Parravicini et al., 2015). Such a niche expansion can be facilitated by adaptive evolutionary processes in the novel distribution area (Reznick et al., 2019;Szűcs et al., 2017), different biotic interactions, or from preadaptation to conditions not available (anymore) in the species' native range but available for the introduced populations (Guisan et al., 2014;Pearman et al., 2008). Hence, predicting future species distributions, by using only the native climatic niche, might severely underestimate the species' geographic potential; this illustrates the importance of evaluating already existing niche shifts, and thus, the need to integrate non-native occurrences when predicting the potential future range.
The exact global distribution of the two mosquitofish species has been difficult to establish for several reasons. First, both species are morphologically very similar and also hybridize, which makes identification challenging and results in many miss-identifications (Scribner & Avise, 1993;Walters & Freeman, 2000). Second, both species were listed as a subspecies of G. affinis until 1988 (Pyke, 2008;Wooten et al., 1988). This complicates the evaluation of historical introduction events and results in a number of erroneous species identifications even today. For example, earlier literature reported the presence of G. affinis in Europe (Innal & Erk'akan, 2006;Krumholz, 1948), while more recent studies could only prove the presence of G. holbrooki (Santi et al., 2020;Sanz et al., 2013;Vidal et al., 2010). This is also the case for genetic studies of mosquitofish in Australia, where only G. holbrooki could be detected so far (Ayres et al., 2010(Ayres et al., , 2013. The situation is less clear for Asia, Latin America, and Africa, where we lack large-scale genetic studies. However, regional molecular surveys suggest only one species, G. affinis, to be common in several southeast Asian countries (e.g., mainland China, Taiwan, Malaysia, Myanmar; Chang et al., 2019;Gao et al., 2017;Kano et al., 2016;Ouyang et al., 2018;Walton et al., 2016).
This global lack of clarity regarding the distribution of both species makes it challenging to evaluate the environmental requirements for both species separately. On the other hand, both species seem to pose a similar threat to indigenous fauna (Pyke, 2005(Pyke, , 2008, have a close taxonomic relationship (Lydeard et al., 1995), and similar ecologies (Walton et al., 2012). Moreover, there is an urgent need to identify the potential distribution range of these highly successful IAS as early as possible and thus establish strategies to detect and prevent further spread. In order to address these needs, our study sought to identify the areas climatically suitable for the establishment and possible spread of mosquitofish, while also considering potential niche shifts of established populations. More specificallyand in order to address different levels of certainty of species identification-we considered different taxonomic levels to address three interrelated questions: 1. In a first approach we treated both species together (i.e., combined species approach) and asked whether invasive mosquitofish conserve their climatic niche between native and introduced ranges. Therefore, we merged occurrence data from G. affinis and G. holbrooki and predicted that-due to an enormous introduction effort at the beginning of the last century-the global invasive range does already cover the entire climatic niche of the native range.
2. In a second, more speculative approach, we adopted a speciesspecific approach, assuming that we can extrapolate evidence from genetic studies to all mosquitofish occurrences in the respective region. We predicted these models to reveal species-specific climatic preferences that lead to different regional establishment probabilities for the two species.
3. Finally, we applied species distribution models to (a) identify areas prone to invasion under current climatic conditions and (b) to project climate change-induced range shifts of both species. Such models will improve the assessments of species' invasive potential and guide future management actions.

| ME THODS
We first examined and identified the climatic conditions under which the species occur (or have been successfully established) and compared these climatic niches between native range and non-native range as well as between G. affinis and G. holbrooki. The climatic niche is defined (here) as the range of climatic conditions under which a species occurs (part of the niche space). Niche overlap indicates the range of climatic conditions under which two species can both occur; niche unfilling and niche expansion indicate the proportion where only one of the two species occurs, or-with regard to invasive species-niche unfilling refers to the climatic conditions under which the species occurs in the native range but not (yet) in the non-native range. In analogy, niche expansion represents the range in which the species occurs in the non-native range but not in the native range (Guisan et al., 2014).
Specifically, we compared (a) the native range niche with the non-native range niche (combined species approach) and (b) the species-specific niches by continent. Based on the niche comparisons using the framework implemented in the R package ecospat (see below), we then built an ecological niche model and projected the global climatic suitability for the two Gambusia species and thus the potential distribution of the species and estimate a future dispersal potential under changed future climate conditions.

| Species distribution data
We obtained information on the global distribution for G. affinis and G. holbrooki from the Global Biodiversity Information Facility (GBIF, 2020), which covers information from several biodiversity databases including fishbase (Froese & Pauly, 2021). We further screened the existing literature to include additional records of Gambusia affinis/G. holbrooki establishment. From this, a total of 64,945 raw records were obtained for both species (Table A1). We used the R package CoordinateCleaner (Zizka et al., 2019) and the spThin package (Aiello-Lammens et al., 2015) to flag potentially erroneous coordinates. Finally, a manual plausibility check was performed to validate the final dataset. We considered only one occurrence record per grid cell (spatial resolution of 10 arcmin as for climatic conditions) even if more than one occurrence was recorded, resulting in 8419 presences (4826 in the native range and 3593 in the non-native ranges-combined species approach; and for the species-specific approaches: for G. affinis 3745 native occurrences and 655 non-native occurrences in North America, 139 occurrences in Asia, and for G. holbrooki 1290 native range occurrences in North America, 894 occurrences in Europe, and 1370 in Australia).
In a preliminary analysis we further cleaned the data based on the year of sampling so that they were in line with the climatic data (i.e., we only considered occurrences between 1970 and 2000).
However, the spatial patterns in the occurrences vary between decades and do not reflect temporal changes in actual distribution but rather temporal changes in sampling effort. For example, almost all occurrences from the Iberian Peninsula are from the period after 2000, but it is well known that G. holbrooki was also widespread there in the period 1970-2000 (Krumholz, 1948;Pyke, 2008;Vidal et al., 2010). The massive and widespread introduction campaigns of Gambusia ended in the 1980s or earlier. Thus, the absence of the species in this area is due to a lack of sampling/reporting to GBIF and could not be explained by recent climatic changes. Accordingly, after careful consideration, we decided to continue working with the complete dataset, arguing that we would introduce a sampling bias into the data that is not justified by the advantage of matching the occurrence data to the related climate data.
Occurrences that were north of the native distribution area but connected by river systems to the south were also classified as native, as we assumed that natural dispersal processes were just as likely as human introductions. Accordingly, the native range of G. holbrooki spans from Alabama, east into Florida, and north along coastal drainages to New Jersey ( Figure A3), whereas G. affinis occurs from Alabama, west into New Mexico, and the Gulf drainages of eastern Mexico ( Figure A4). A hybrid zone between the two closely related species can be found in the Mobile Bay, Alabama region (Wilk & Horth, 2016). In this area of overlapping distribution, we adopted the species information from the raw data (see Figures A1-A9 for details).

| Environmental variable selection
We used 19 climatic variables at 10-arcminute resolution from the WorldClim database (Fick & Hijmans, 2017). To deal with strong collinearity among climatic variables, we calculated Pearson's correlation coefficients (r) between the 19 climatic variables and constructed a cluster dendrogram ( Figure A10). Based on a threshold of r ≥ |0.8| (Elith et al., 2006;Franke, 2010;Mateo et al., 2013) nine groups of intercorrelated variables were found. From five of these groups, we have chosen one representative that we consider to be ecologically most meaningful for the distribution of the fish species (Gao et al., 2017;Riesch et al., 2018;Santi et al., 2020). The climatic variables that we selected for the analyses were temperature seasonality (bio4), maximum temperature of the warmest month (bio5), minimum temperature of the coldest month (bio6), annual precipitation (bio12), and precipitation seasonality (bio15). We did not include additional precipitation variables (i.e., precipitation in the coldest [bio19] and the driest periods [bio14, bio17]), as well as the mean diurnal range (bio2) and isothermality (bio3), because we assume that diurnal differences in air temperature were thermally buffered by the water and, hence, less relevant for aquatic species.

| Comparison of the climatic niches
To quantify potential shifts in the niches of G. affinis and G. holbrooki, we used the Centroid, Overlap, Unfilling, Expansion (COUE) framework of Guisan et al. (2014) to decompose niche changes into centroid shifts, degree of overlap, and amounts of unfilling and expansion. We used the R package ecospat (Broennimann et al., 2020;Di Cola et al., 2017) to investigate the species' distribution in the niche space. The resulting smoothed occurrence densities were plotted into the ordination space of PCA (based on five bioclimatic variables: bio04, bio05, bio06, bio12, bio15) to visualize the position of within environmental space (i.e., the realized/occupied climatic niche). Specifically, we compared native range and the non-native range niches for both species together (combined species approach) and for each species separately (species-specific approach). For each comparison of two niches, we calculated niche unfilling, defined as the percentage of the first niche covered by the second niche. In our case, a high unfilling means that a large part of the native range niche is unoccupied by the species in the non-native range, or, when comparing two species, that there is a large range of environmental conditions under which the first species occurs but not the second.
In addition, we calculated niche expansion, defined as the percentage of the second niche covered by the first niche. For our data, a high niche expansion means that the species occurs in the nonnative range under new conditions under which it does not occur in the native range, or, when comparing two species, that there is a large range of environmental conditions under which the second species occurs but not the first. Absolute overlap between the two given niches was further calculated based on the position of the occurrence densities using the D metric of Schoener (1968), which ranges between 0 (no overlap) and 1 (complete overlap).

| Species distribution modeling
We projected the global habitat suitability for G. affinis and G. holbrooki (combined species approach) under current and future climatic conditions based on the ecological niche modeling approach.
The ecological niche modeling was performed with an ensemble forecasting approach incorporating six state-of-the-art niche modeling algorithms (ANN-artificial neuronal networks, GAMgeneralized additive models, GBM-generalized boosted models, GLM-generalized linear models, FDA-flexible discriminant analysis, and RF-random forest) and executed in the R environment (R Development Core Team, 2019) using the biomod2 package (Thuiller et al., 2020). In the ensemble modeling approach, single model results are merged into a consensus model (here: weighted average), which is then considered to be a more robust estimator as it reduces uncertainties due to the choice of algorithm (Araújo & New, 2007).
In the ensemble model, we consider all presence-absence algorithms available in biomod2 (i.e., ANN, GAM, GBM, GLM, FDA, RF) and excluded the presence-only algorithms (SRE) and presence background algorithms (Maxent), as we believe that due to the intensive introduction history, missing distribution data have a high information value and can be evaluated as absence data.
Ten thousand pseudo-absences were chosen at random but excluding the area close to observed presences of the species in order to avoid pseudo-replication, as close points tend to show similar environmental conditions/same niche (disk strategy implemented in the biomod2 package). The models were run using the following single algorithm parameters: Running the artificial neuronal networks (ANN) we used five cross-validation to find the best size and decay parameters, and Continuous modeling results were transformed into binary results, using the equal sensitivity and specificity threshold (Liu et al., 2005). We subtracted the binary maps for the current conditions from the projected future distribution (2081-2100) to quantitate the impacts of climate change on range size changes. Based on these binary results (climatically suited or unsuited) under current and future climatic conditions, we thus identified areas that were (a) currently suited but projected to become climatically unsuitable in future (po- We also looked at the variables' contribution to the ensemble model of the considered five climate variables to identify important driving factors in the invasion process. In an analogous way, we have used the species-specific occurrence data to model the speciesspecific climatic suitability ( Figures A18 and A19). Based on this (more speculative) approach, we identified areas that are climatically suited (according to our model) for both or only one of the two species, respectively ( Figure A21).

| Projections of potential distributions
The current distribution of both mosquitofish species (combined Only in South America, a southern range expansion is to be expected. The greatest potential for range expansion can be expected in Europe. Regardless of the climate scenario considered, habitats in central Europe and southern United Kingdom will represent suitable habitats in the future. Assuming a moderate (SSP 3.70) or extreme (SSP 5.85) climate change scenario, the suitable habitat could even increase as far as southern Scandinavia and western Russia by the end of the century. In contrast, northern distributions in the F I G U R E 1 (a) Combined analysis of Gambusia affinis and G. holbrooki niche shift (combined species approach) in global environmental space, derived from principal component analysis on the climate predictors. Solid contour line represents available climates and dashed line the 50% most frequent available climate. Blue shaded area represents the niche area occupied in both native and invasive range; green area is the unfilled niche in the invasive range (relative to the native range), and pink shows the expansion area. The gray shading within these contours (black) correspond to the densities of occurrence records within the occupied climatic space of the latter niche (here: non-native range niche). (b) Correlation circle indicates the weight of the selected climatic variables on the niche space as defined by the first two principal component axes (explaining 79.64% of the variance in the set of five predictor variables); bio4 = temperature seasonality, bio5 = maximum temperature of the warmest month, bio6 = minimum temperature of the coldest month, bio12 = annual precipitation, bio15 = precipitation seasonality TA B L E 1 Pairwise niche overlap indices (Schoener's D) of Gambusia affinis and G. holbrooki between native and invasive ranges

F I G U R E 2
Projection of the realized niches in climatic space (species-specific approach), comparing (a) populations within the native range of Gambusia holbrooki and G. affinis; native and invasive range of G. affinis in (b) North America and (c) Asia; native and invasive range of G. holbrooki in (d) Europe and (e) Australia, as well as (f) the two invasive ranges in Europe and Australia. Solid contour line represents available climates and dashed line the 50% most frequent available climate. Blue areas symbolize niche overlap; green area is the niche exclusively filled by the first-mentioned species (i.e., "unfilled" by the second-mentioned species), and pink shows the "expansion" area, that is, the climatic niche space solely occupied by the second-mentioned species. The gray shading shows the smoothed occurrence density of the latter mentioned niche  Table A2).

| DISCUSS ION
In the first part of our paper, we explored niche dynamics of the highly invasive G. affinis and G. holbrooki in a combined species approach. In accordance with our associated prediction, we found that invasive mosquitofish occupy geographic areas that share the full range of climatic conditions occupied in the native range (only 1% unfilling). Moreover, our results suggest a slight expansion of the climatic niche by 8% during, or subsequent to, invasion of both species.
This niche expansion becomes even more evident in our second, species-specific approach: When comparing invasive populations of G. holbrooki in Europe with their native distribution range, we found a pronounced niche shift with a large degree of both expansion (81%) and low overlap with the native niche. Using a species-specific approach for European populations is reasonable, because we have sufficient evidence from genetic studies that only G. holbrooki is present (Santi et al., 2020;Sanz et al., 2013;Vidal et al., 2010). The genetic evidence for the exclusive introduction of G. holbrooki is also good in Australia (Ayres et al., 2010(Ayres et al., , 2013, and here we also find a niche expansion, albeit less pronounced, compared to Europe (47% niche expansion). However, we found the most pronounced niche  (Havel et al., 2015;Pyšek et al., 2008). In general, the observed expansion, especially in Asia, is not surprising due to the outstanding ability of both species to tolerate a wide range of environmental conditions, being extremely flexible in terms of their habitats (Pyke, 2008), diet (Pirroni et al., 2021) and abiotic conditions (Cherry et al., 1976;Odum & Caldwell, 1955). Mosquitofish respond to multiple interacting environmental factors by seemingly adaptive life-history shifts Santi et al., 2020), which promote invasiveness and facilitate the colonization of new environments (Hendry, 2016).
The observed niche expansion suggests that their realized niches in North America actually do not encompass their entire physiological and ecological ranges (Rosenfield, 2002). Indeed, native species' distribution is often limited by biotic constraints (e.g., predation, competition, parasitism) and/or by biogeographical barriers (Moore et al., 2007;Zaret, 1980). While the native distribution range of G. holbrooki is limited primarily by the Atlantic Ocean, the Gulf of Mexico, and/or competition from G. affinis, G.
affinis encounters tremendous competition from other Poeciliidae at the southern range boundary, as Mexico already harbors >100 species of Poeciliidae, including 25 species of Gambusia (GBIF, 2021).
If these constraints were removed, both mosquitofish species seem capable of occupying a much wider geographical and ecological range of habitats.
We detected these niche shifts even though we used a very conservative approach when defining the native range (i.e., we classified northern populations in eastern North America as native). Niche shifts following biological invasions have been repeatedly described (Broennimann et al., 2007;Gallagher et al., 2010), although a recent meta-analysis for 434 invasive plant and animal species generally demonstrated that there is very limited niche expansion between native and introduced ranges (Liu et al., 2020). However, the occurrence of niche expansion was most evident in aquatic species, suggesting they might be more capable of invading diverse environments and highlighting the importance of aquatic habitats for conservation and species management (Liu et al., 2020). In the case of European mosquitofish, our study indicates that the niche shift occurred despite their low genetic diversity in Europe (Santi et al., 2020) and a short time span since the introduction at the beginning of the 20th century (Krumholz, 1948). We, therefore, can assume that the observed climatic niche shift of both mosquitofish might be a combination of microevolutionary change and/or adaptive plasticity, that is, shifts in the realized climatic niche within the broad fundamental niche of G. affinis and G. holbrooki. The latter matches recent findings regarding other aspects of their phenotype (Santi et al., 2020) and would represent a preadaptation to conditions that are not present in the native area (Cadotte et al., 2018;Petitpierre et al., 2012), indicating an intrinsic capacity to be successful invaders of novel environments.
Based on our finding that both mosquitofish expanded their realized niche globally beyond their native niche(s), we have described the potential range of G. affinis and G. holbrooki using the pooled occurrence data from native and invasive range. This represents a first step toward a risk analysis of both species, which are already recognized as highly problematic global invaders (Lowe et al., 2000). Our projections under future climatic scenarios are alarming and show that (particularly) large areas of Central Europe are predicted to become climatically suitable for mosquitofish in the future. So far, there have been no global models for the probability of establishment, and so it is not surprising that the risk of invasion is assessed differently by different European countries. For example, a recent horizon scan for potential future IAS threatening Great Britain biodiversity did not consider mosquitofish (Roy et al., 2014), while they are on the watch list in Germany, but without specific actions or plans for management (Nehring et al., 2010). The main reason for the severe risk of future mosquitofish establishment are elevated winter temperatureswhich so far have prevented the establishment of G. holbrooki in central Europe (Kinzelbach & Krupp, 1982). This highlights the urgent need to devise appropriate species management plans for the areas and countries predicted to be affected in the near future.
A general shortcoming of species distribution models of aquatic species is that they often refer to terrestrial climate scenarios.
Unfortunately, the resolution of the aquatic predictors (Domisch et al., 2015) is not yet sufficient to properly enable the analyses we conducted here. Gambusia holbrooki and G. affinis are species that occur in the smallest puddles and roadside ditches, about which no robust global environmental information is yet available. However, these small habitats are also strongly influenced by the surrounding air temperatures due to their low water masses/volumes, which is why we are convinced that the terrestrial climate data used here are suitable predictors.
In conclusion, our results show that a combination of niche conservatism and niche expansion facilitate a large climatic tolerance that helps explain the observed invasive success of both species in several parts of the world and indicate that there is potential for further range expansion in the face of global warming. The control and eradication of mosquitofish is often promoted by both governmental and scientific authorities (Fryxell et al., in press;Pyke, 2005Pyke, , 2008. However, limited evidence exists on feasibility of mosquitofish removal from natural environments in which they are established but not native (Brookhouse & Coughran, 2010;Cano-Rocabayera et al., 2019). Contingency plans should, therefore, focus on prevention, especially in regions with suitable climatic conditions. One important tool is to reduce anthropogenic pressures (e.g., habitat modification, dam construction, pollution loads) and put effort into the restoration of disturbed habitats, because it is becoming increasingly evident that IAS flourish primarily in anthropogenically disturbed river systems (Lee et al., 2017). (Table A3, Figures A11, A12, A13, A14, A15, A16, A17.)

ACK N OWLED G M ENTS
We thank Christian Hof for a very useful discussion of ideas at the beginning of the study. We also thank Zhixin Zhang and anonymous reviewers for valuable comments and suggestions on earlier drafts of the manuscript.

CO N FLI C T O F I NTE R E S T
The authors have not declared any conflict of interest.

F I G U R E A 2
Occurrence records per grid cell for the combined invasive range of Gambusia affinis and G. holbrooki (Combined species approach; 6889 data points)

F I G U R E A 3
Occurrence records per grid cell for the native range of Gambusia holbrooki (3279 data points)

F I G U R E A 4
Occurrence records per grid cell for the native range of Gambusia affinis (8593 data points)

F I G U R E A 5
Occurrence records per grid cell for the invasive range of Gambusia affinis in North America (1167 data points)

F I G U R E A 6
Occurrence records per grid cell for the invasive range of Gambusia holbrooki in Europe (1488 data points)

F I G U R E A 7
Occurrence records per grid cell for the invasive range of Gambusia holbrooki in Australia (3053 data points)

F I G U R E A 8
Occurrence records per grid cell for the invasive range of Gambusia affinis in southeast Asia (185 data points)

F I G U R E A 9
Global occurrence records per grid cell for Gambusia affinis and G. holbrooki (18,508 data points)

A PPE N D I X 2
Environmental variable selection F I G U R E A 1 0 Cluster dendrogram with the Pearson's correlation coefficients (r) between the 19 climatic variables provided by worldclim. Based on a threshold of r ≥ |0.8 | (blue line) nine groups of intercorrelated variables were found. From five of these two groups we have chosen one representative that we consider to be ecologically most meaningful for the distribution of the fish species (blue arrows). We did not include additional precipitation variables (i.e., precipitation in the coldest [bio18, bio19] and the driest periods [bio14, bio17]), as well as the mean diurnal range (bio2) and isothermality (bio3), because we assume that diurnal differences in air temperature were thermally buffered by the water and, hence, less relevant for aquatic species Variable contribution to the combined species ensemble model

A PPE N D I X 7
Country-specific invasion risk index for G. affinis/G. holbrooki F I G U R E A 2 0 Variable importance of different models in the combined species approach TA B L E A 2 Categorized establishment risk that at least one of the two species (Gambusia affinis/G. holbrooki) finds suitable climatic conditions in large parts of the respective country (under current and predicted climatic conditions of the 2081-2100 period; four SSPs). Relative climatic habitat suitability was calculated for each pixel in a country and averaged it over the total area of the country: minimal risk (0-14); low risk (0.15-0.24), moderate risk (0.25-0.49), and high risk (0.50-1.0). For large area countries that are climatically heterogeneous (e.g., China), this approach may result in low probabilities, although the probability of local establishment is very high. These countries are marked with an asterisk (*)