The Gulf toadfish goes south: assessing the potential habitat suitability of Opsanus beta in response to climate change

. Climate change can directly influence invasive species range, and may reflect in greater habitat suitability, consequently increasing losses in native biodiversity. Habitat suitability modelling can help identify priority areas for conservation and rapid response for invasive taxa. This study aimed to model the habitat suitability of Opsanus beta , an invasive fish species, in South America under current and future climate change scenarios. Three modelling techniques, Generalized Linear Model (GLM), Maximum Entropy Algorithm (Maxent), and Random Forest (RF), were applied, and their results were compared and then put together in an Ensemble model to visualize habitat suitability. We found that the current habitat suitability of O. beta in South America is relatively low, with its occurrences mostly related to coastal port areas, suggesting that the species was introduced from ballast water. The increase in mean temperature at mean depth and temperature range were the most relevant variables influencing the species habitat suitability. The gradual increase in habitat suitability for the 2100 climate change scenarios, particularly in more severe climate change scenarios, such as the RCP 8.5, was also observed. The study highlights the critical need to use habitat suitability models to mitigate the impact of future climate change scenarios on invasive species. The findings call for preventive measures to be taken in time to prepare areas that may become the target of dispersal and establishment of O. beta .


Introduction
An urgent threat to the world today is posed by invasive species, which disrupt ecological balance by displacing and even driving native species to extinction.This threat has become even more evident in the context of climate change, which directly influences the potential habitat suitability for these species in various regions (Holcombe, Stohlgren, & Jarnevich, 2010).Moreover, the number of invasive species that are relocated due to human activities is rapidly increasing, making them one of the most dangerous generators of changes and losses in native biodiversity.Therefore, it is essential to better understand the potential invasions and expansions of these invasive species (Seebens, Schwartz, Schupp, & Blasius, 2016).
Shipping trade has proven to be a powerful route for the introduction of exotic species (Wan, Shi, Nie, Chen, & Wang, 2021) and predicting the proliferation of such species that have arrived through this type of route will enable a more effective response and proper management.Habitat suitability modelling (HSM) can identify the most potentially vulnerable areas, as rapid detection and response are requirements for the management of invasive species (Crall et al., 2013;Goldsmit et al., 2018), particularly those that have not yet occupied their potential geographical range.The HSM can provide the potential ecological niche of these species based on environmental variables in the region where the species is considered native and wellestablished, with predictions on the invaded regions (Srivastava, Lafond, & Griess, 2019).Furthermore, we need to consider the impact of climate change on these potentially suitable areas, which may even favor the establishment of certain invasive species (Chai, Zhang, Nixon, & Nielsen, 2016).
One species that has become recurrent on the South America coast, especially on the Brazilian coast is Opsanus beta (Goode & Bean, 1880), with reports in the Santos Estuary, where the largest port in Latin America is located, and with dispersal being recorded for other areas through citizen science (Carvalho et al., 2022).Mature individuals have been sampled, indicating reproductive activity in the region (Tomás, Tutui, Fagundes, & Souza, 2012).There are also reports along the coasts of Rio de Janeiro and Paraná, Brazilian states that have estuarine and port areas, suggesting that its introduction may have come from ballast water (Caires, Pichler, Spach, & Ignácio, 2007;Cordeiro et al., 2020), although there is the possibility of introduction by other human-aid displacement (Carvalho, Ferreira Junior, Fávaro, Artoni, & Vitule, 2020), such as the aquarium trade (Carvalho et al., 2022).This species is native to the Western Central Atlantic and is known for its occurrence in the Gulf of Mexico.It is classified as a resilient, territorial, and aggressive species that feeds on fish, crustaceans, and mollusks, with a preference for warm waters, what emphasizes the need to understand its potential distribution and identify the areas at risk of proliferation (Andrade-Tubino, Salgado, Uehara, Utsonomia, & Araújo, 2021).In light of recent warnings given by Miller, Carvalho, Figueiredo, Aguilar-Perera, and Vitule (2023), it is imperative to vigilantly track the geographical dispersion of the species in regions where it is not native and to investigate potential shifts in the species distribution precipitated by future climate change scenarios.Hence, this study aims to understand the habitat suitability of O. beta on the South American coast, using environmental data from the region where the species is considered native, to identify the areas with the highest potential for the occurrence of the invasive species, not only in the present but also under different scenarios of climate change, from mild to severe, with the aim of assisting in the proper management and control of this invasion while there is still time.

Material and methods
The region delimited for the present study was South America, specifically focusing on Brazilian coastal regions where O. beta has been recorded as an invasive species.Data to model the species distribution was collected from various sources including capture, human observations, literature publications, through the Global Biodiversity Information Facility online database (GBIF.org),with the rgbif R package (Chamberlain, 2019), being selected a total of 5000 records of the species in its native area.The data was cleaned to eliminate duplicates or records lacking longitude and latitude information and converted to degrees using the WGS84 global reference system.Environmental variables of the marine environment were obtained from the Bio-Oracle global database (Tyberghein et al., 2012;Assis, Tyberghein, Bosch, Verbruggen, Serrão, & De Clerck, 2018) through the R package sdmpredictors (Bosch, 2018), with 5 arcmin resolution and selected based on availability for both current conditions and Intergovernmental Panel on Climate Change (IPCC) Representative Concentration Pathway (RCP) scenarios.The following variables were considered: temperature mean at mean depth, temperature range at mean depth, temperature mean at surface, temperature range at surface, salinity mean at surface, salinity range at mean depth, and current velocity mean at surface.The correlation between the variables was tested using Pearson's correlation, retaining only those with a correlation value lower than 0.7.
The models were fitted using records from areas that O. beta is considered native where the species is in balance with its environment (Figure 1), and predictions were made for South America coast, evaluating habitat suitability for the current time and three scenarios in the period from the years 2090 to 2100: mild (RCP 2.6), intermediate (RCP 6.0), and severe (RCP 8.5).Three modelling techniques were employed: Generalized Linear Model (GLM), Maximum Entropy Algorithm (Maxent), and Random Forest (RF), using the R sdm package (Naimi & Araujo, 2016).5000 pseudo-absence records were randomly generated within the species recorded range and the models were tested using three bootstraps of the data for each model technique.The accuracy of the model was evaluated using the Area Under the Receiver Operating Characteristic Curve (AUC) and True Skill Statistics (TSS), and the performance of each model was visualized through a Receiver Operating Characteristic Curve (ROC curve).Overall importance of every variable was measured and the average response curves of each modelling technique for the selected environmental variables was visualized.
Ensemble models combining GLM, Maxent and RF results were generated by averaging the predictions, balanced by True Skill Statistics (TSS) values to visualize suitability in a more reliable way, generating maps of habitat suitability for the current time and three RCP scenarios in 2100, also considering the difference in suitability between RCP 8.5 and the current time, that was calculated and visualized as a map.To characterize ecological niche considering only environmental variables, the habitat suitability of the species was mapped into a two-dimensional environmental space considering temperature mean, temperature range and salinity mean.All analyzes were performed using the R software (R Core Team, 2022).(Chamberlain, 2019).

Results
The correlation between all environmental variables was examined, as illustrated in Figure 2. At the time all models for predicting the habitat suitability of O. beta were successfully executed.The results for each modelling technique in terms of AUC and TSS values are presented in Table 1.The Random Forest method emerged as the top performer, boasting the highest TSS value and lowest deviance.While the results from various models were relatively similar, with Maxent and Random Forest exhibiting the highest AUC, with the former presenting a higher deviance.The predictions were made considering all the models and gathered in an ensemble model, which was used to visualize the habitat suitability in the present time and different RCP's scenarios, with an increase of the suitability in more severe scenarios of climate changes, as in the RCP 8.5 (Figure 3).To clarify the distinction between the current suitability and the worst-case scenario, the difference in suitability was analyzed (Figure 4), revealing a substantial change in habitat suitability.The response curves of O. beta habitat suitability were explored by model technique (Figure 5), with all models showing a high relationship of the species with the increase of mean temperature in a mean ocean depth (BO2_tempmean_bdmean).For Random Forest the same happened with temperature range at a mean ocean depth (BO2_temprange_bdmean), while other variables did not show an evident relationship with the species.Seeing that, the importance of the relative variable of each variable was compared, considering all techniques altogether (Figure 6).The ecological niche of the species mapped in two-dimensional environmental space showed that O. beta habitat suitability is higher when the temperature values increase independently of the temperature range.The same happens when temperature is mapped with salinity (Figure 7), although there is a predominance of greater habitat suitability in high salinity values.Low salinity measurements did not stop the increase in suitability.(Chamberlain, 2019) and environmental data retrieved from the Bio-Oracle global database (Tyberghein et al., 2012).(Chamberlain, 2019) and environmental data retrieved from the Bio-Oracle global database (Tyberghein et al., 2012).
Figure 7. Two-dimensional representation of the Opsanus beta ecological niche only considering environmental variables, respectively the relationship of temperature mean at mean depth with temperature range at sea surface (A), and the relationship of temperature mean at mean depth with salinity mean at sea surface (B).The higher the value of the niche, the greater is the habitat suitability.
Occurrence data retrieved from rgbif R package (Chamberlain, 2019) and environmental data retrieved from the Bio-Oracle global database (Tyberghein et al., 2012).

Discussion
We have successfully modelled the habitat suitability of O. beta by applying three modelling techniques, Generalized Linear Model (GLM), Maximum Entropy Algorithm (Maxent), and Random Forest (RF).The results of all three approaches were comparable, the RF model exhibited the highest True Skill Statistic (TSS) value, which is used to measure a model's ability to differentiate between occurrences and absences (Somodi, Lepesi, & Botta-Dukát, 2017), and the lowest deviance, while both RF and Maxent models showed relatively high Area Under the Curve (AUC) values, RF still emerged as the top performer, since Maxent presented a high deviance value, and GLM expressed a lower TSS value and a slightly greater deviance than RF.
It is worth noting that the current habitat suitability of O. beta in South America is relatively low, suggesting that the species was probably introduced from ballast water, as previously reported (Cordeiro et al., 2020).Consequently, its occurrences are mostly related to coastal port areas, rather than regions with greater suitability.Nevertheless, the present areas with greater suitability are those close to the coast in estuarine regions, particularly in southeastern Brazil, as demonstrated in Figure 3. Opsanus beta has a high probability of being an established species in estuarine regions, subjecting other estuarine regions to the threat of its invasion through the continuity of its dispersion by ship traffic (Tomás et al., 2012).Furthermore, the species is well stabilized in the Santos Port, one of the largest recipients of ships in Latin America, and the success of its establishment in these estuaries is related to the warm water surrounding the mangroves and the absence of predators, such as some species of groupers and marine mammals (Andrade-Tubino et al., 2021), with the species even nesting in human litter (Barimo, Serafy, Frezza, & Walsh, 2007), that are abundant in port regions.
The response curves of O. beta habitat suitability in Figure 5 demonstrated a high relationship with the increase in mean temperature at mean depth for all models.However, for the RF model, the same relationship was observed with temperature range, with these variables being the most relevant ones, as shown in Figure 6.Interestingly, in the two-dimensional representation of the ecological niche of the species, the habitat suitability was higher when temperature values increased, regardless of the temperature range, and the same pattern was observed when temperature was mapped with salinity.This may be related to the fact that the species thrives in high temperatures, even in places with low and high temperature variations (Figure 7).Similarly, the invasion of the species is indeed disturbing, not only because it adapts very well to high temperature variations, but also to a wide range of salinity values.The preference for high temperature makes sense, as they naturally occur in warmer embayment in Florida and along the Gulf of Mexico coast (Malca, Barimo, Serafy, & Walsh, 2009), and this could explain the fact that the species has established itself in port regions shortly after the invasion.
However, the most concerning finding of our study is the gradual increase in habitat suitability for the 2100 climate change scenarios, particularly in scenarios of more severe climate changes.The difference between the visualization of the current and the RCP 8.5 habitat suitability ensemble models, as illustrated in Figure 4, revealed a substantial change in habitat suitability.This indicates that if a severe scenario such as RCP 8.5 occurs, there is a stronger possibility of the species dispersal to other coastal areas, including the northeast of South America, an area where the species has not yet been observed.In addition to the likely increase in temperature in the coming years facilitating the occurrence of the species, it is important to consider that the impact of its possible proliferation could be costly on native species, since O. beta is a territorial species, has high parental care, and a faster growth than other species of the genus (Malca et al., 2009), also presenting a great potential for invasion due to being an opportunistic and predatory species (Miiller, Carvalho, Figueiredo, Aguilar-Perera, & Vitule, 2023).
These results highlight the critical need to use habitat suitability models to mitigate the impact of future climate change scenarios on invasive species, which is increasing in biodiversity assessments (Araújo et al., 2019).This will enable preventive measures to be taken in time, acting as a strategy to prepare areas that may become the target of dispersal, and establishment of these species (Seebens et al., 2016).Without a doubt, given the lack of future predictions of some variables, future research must explore different environmental variables at distinct spatial scales even for present times.

Conclusion
Our study effectively modeled O. beta habitat suitability, identifying the Random Forest model as best performer.We found that this species thrives in high temperatures and varied salinities, subjecting other estuarine regions to the threat of its invasion.The rising suitability under 2100 climate scenarios suggests potential spread along the South American coast.These insights underscore the importance of habitat models in forecasting invasive species threats under climate change and emphasize refining environmental variables at diverse spatial scales.Our research aids in devising strategies to mitigate invasive species impacts and safeguard coastal ecosystems and dependent communities.

Figure 1 .
Figure 1.Distribution of the Global Biodiversity Information Facility records on Opsanus beta (blue points), with the native range species used to fit the model and the predicted range for habitat suitability.Data retrieved from rgbif R package (Chamberlain, 2019).

Figure 2 .
Figure 2. Pearson correlation for all environmental variables collected from Bio-Oracle and used to model the habitat suitability of Opsanus beta.Data retrieved from the Bio-Oracle global database (Tyberghein et al., 2012).

Figure 3 .
Figure 3. Ensemble model visualization of Opsanus beta habitat suitability in the present time and three RCP's scenarios from IPCC, with suitability values ranging from 0 to 1. Occurrence data retrieved from rgbif R package(Chamberlain, 2019) and environmental data retrieved from the Bio-Oracle global database(Tyberghein et al., 2012).

Figure 4 .
Figure 4. Difference between the current and RCP 8.5 habitat suitability ensemble models visualization of Opsanus beta in the South America coast, the higher values represent a habitat suitability increase.Occurrence data retrieved from rgbif R package(Chamberlain, 2019) and environmental data retrieved from the Bio-Oracle global database(Tyberghein et al., 2012).

Figure 5 .Figure 6 .
Figure 5. Response curves of Opsanus beta habitat suitability for each modelling technique, respectively Generalized Linear Model (A), Maximum Entropy Algorithm (B) and Random Forest (C).BO2_curvel_mean_ss = Current velocity mean at surface; BO2_salinitymean_ss = Salinity mean at sea surface; BO2_salinityrange_bdmean = Salinity range at mean depth; BO2_salinityrange_ss = Salinity range at sea surface; BO2_tempmean_bdmean = Temperature mean at mean depth; BO2_tempmean_ss = Temperature mean at sea surface; BO2_temprange_bdmean = Temperature range at mean depth; BO2_temprange_ss = Temperature range at sea surface.Occurrence data retrieved from rgbif R package (Chamberlain, 2019) and environmental data retrieved from the Bio-Oracle global database (Tyberghein et al., 2012).

Table 1 .
Model performance values for the three approaches applied to model the habitat suitability of Opsanus beta.GLM = Generalized Linear Model; Maxent = Maximum Entropy Algorithm; RF = Random Forest; AUC = Area Under the Receiver Operating Characteristic Curve; TSS = True Skill Statistics.