Warming Threatens to Propel the Expansion of the Exotic Seagrass Halophila stipulacea

The spread of exotic species to new areas can be magnified when favored by future climatic conditions. Forecasting future ranges using species distribution models (SDMs) could be improved by considering physiological thresholds, because models solely based on occurrence data cannot account for plasticity due to acclimation of individuals to local conditions over their life-time or to adaptation due to selection within local populations. This is particularly relevant for the exotic seagrass Halophila stipulacea, which colonized the Mediterranean Sea a century ago and shifted its thermal niche, coping with a colder regime. Here, we used two hybrid models combining correlative SDMs with the thermal limits for growth of native and exotic H. stipulacea populations to predict the distribution of the species in its native (Indian Ocean and Red Sea) and exotic ranges (Mediterranean Sea and Caribbean Sea) under two scenarios forecasting limited (RCP 2.6) and severe (RCP 8.5) future climate changes by 2050 and 2100. Then, we assessed the differences between hybrid models based on native Red Sea thermal limits (niche conservatism: 17–36°C) and on exotic Mediterranean thermal limits (local adaptation: 14–36°C). At the Mediterranean exotic range, the local adaptation hybrid model accurately agreed with the present distribution of the species while the niche conservatism-based hybrid model failed to predict 87% of the current occurrences of the species. By contrast, both hybrid models predicted similar species distributions for the native range and exotic Caribbean range at present and projected that H. stipulacea will maintain its current worldwide under all future greenhouse gas emission scenarios. The hybrid model based on Mediterranean thermal limits projected the expansion of H. stipulacea through the western Mediterranean basin (except the gulf of Leon) under the most severe scenario (RCP 8.5) by 2100, increasing its distribution by 50% in the Mediterranean. The future expansion of H. stipulacea is related to its capacity to cope with warm waters and it may become a relevant species in the future, particularly under the projected decline of native Mediterranean seagrasses, resulting in important shifts in seagrass communities and overall ecosystem functions.

The spread of exotic species to new areas can be magnified when favored by future climatic conditions. Forecasting future ranges using species distribution models (SDMs) could be improved by considering physiological thresholds, because models solely based on occurrence data cannot account for plasticity due to acclimation of individuals to local conditions over their life-time or to adaptation due to selection within local populations. This is particularly relevant for the exotic seagrass Halophila stipulacea, which colonized the Mediterranean Sea a century ago and shifted its thermal niche, coping with a colder regime. Here, we used two hybrid models combining correlative SDMs with the thermal limits for growth of native and exotic H. stipulacea populations to predict the distribution of the species in its native (Indian Ocean and Red Sea) and exotic ranges (Mediterranean Sea and Caribbean Sea) under two scenarios forecasting limited (RCP 2.6) and severe (RCP 8.5) future climate changes by 2050 and 2100. Then, we assessed the differences between hybrid models based on native Red Sea thermal limits (niche conservatism: 17-36 • C) and on exotic Mediterranean thermal limits (local adaptation: 14-36 • C). At the Mediterranean exotic range, the local adaptation hybrid model accurately agreed with the present distribution of the species while the niche conservatism-based hybrid model failed to predict 87% of the current occurrences of the species. By contrast, both hybrid models predicted similar species distributions for the native range and exotic Caribbean range at present and projected that H. stipulacea will maintain its current worldwide under all future greenhouse gas emission scenarios. The hybrid model based on Mediterranean thermal limits projected the expansion of H. stipulacea through the western Mediterranean basin (except the gulf of Leon) under the most severe scenario (RCP 8.5) by 2100, increasing its distribution by 50% in the Mediterranean. The future expansion of H. stipulacea is related to its capacity to cope with warm waters and it may become a relevant species in the future, particularly under the projected decline of native Mediterranean seagrasses, resulting in important shifts in seagrass communities and overall ecosystem functions.

INTRODUCTION
Temperature regimes, a determinant factor controlling the biogeographic distribution of species (Brown et al., 2004;Dell et al., 2011), are expected to continue to increase due to the rise in anthropogenic emissions of greenhouse gases (Intergovernmental Panel on Climate Change [IPCC], 2014), driving a global redistribution of species with economic and ecological consequences (Pecl et al., 2017). Changes in species phenology and distribution of climatic niches have been documented widely, both on land (e.g., Chen et al., 2011;Gottfried et al., 2012) and ocean (e.g., Sunday et al., 2015;Poloczanska et al., 2016;Wernberg et al., 2016), where marine range shifts tend to be greater and faster than terrestrial ones (Parmesan and Yohe, 2003;Sorte et al., 2010), because of high connectivity and continuity of the marine domain (Burrows et al., 2011).
Species distributions are also altered by human activities that facilitate the spread (Mack et al., 2000;Seebens et al., 2017) and impacts (Geraldi et al., 2020) of exotic species beyond their native biogeographical range, which are expected to increase in the future (Sardain et al., 2019). Climate change is likely to aggravate the degradation of ecosystem services and biodiversity by exotic species (Bennett et al., 2021), which are already substantial (Vilà et al., 2011;Maggi et al., 2015;Doherty et al., 2016;Gallardo et al., 2016;Anton et al., 2019). Similar climatic conditions between native and recipient locations are key for successful biological invasion and establishment (Duncan et al., 2003;Thuiller et al., 2005;Bennett et al., 2021) and, therefore, it is expected that warming might accelerate the spread of exotic tropical species into warm temperate areas and/or increase their abundance, exacerbating their impacts (Bates et al., 2013;Poloczanska et al., 2013;Sorte et al., 2013;Bennett et al., 2021). Besides, in marine systems the impact of exotic species on recipient ecosystems is maximized where the median temperature is slightly (2.2 • C) cooler than the optimal average temperature at the original range of the species (Bennett et al., 2021).
Projections of exotic species distributions under climate change scenarios are critical to anticipate future retractions or expansions in their distributional ranges to guide management actions (Bates et al., 2013;Guisan et al., 2014). Potential future geographical ranges of exotic organisms have been forecasted with species distribution models (SDMs) (Chefaoui and Varela-Álvarez, 2018;Chefaoui et al., 2019;Wilson et al., 2019;Goldsmit et al., 2020;Poursanidis et al., 2020). Correlative SDMs (Guisan and Thuiller, 2005;Elith and Leathwick, 2009) are based on the species' realized distribution, assuming that all relevant ecological processes are captured by the correlative analyses (Guisan and Thuiller, 2005;Buckley et al., 2010;Peterson and Soberón, 2012). However, SDMs do not take into account the mechanisms controlling species distributions and therefore have limitations when making projections for conditions not recorded before, such as new regions or future climates (Elith and Leathwick, 2009;Buckley et al., 2010), which are particularly relevant for the case of exotic species. As an alternative, the physiological tolerance of species has been used to integrate the species' fundamental niche into correlative SDMs to deliver more robust projections (Kearney and Porter, 2009;Buckley et al., 2011). Thus, hybrid of correlative and mechanistic models are increasingly being used to improve marine species distribution projections by incorporating thermal thresholds for survival (e.g., Martínez et al., 2015), growth (Gamliel et al., 2020) or the reproductive phenology of invasive species (Chefaoui et al., 2019).
As temperatures are expected to continue to rise, local populations can also adapt to the changing climate (Matesanz et al., 2010;Nicotra et al., 2010;Hoffmann and Sgrò, 2011;Schlüter et al., 2014;Padfield et al., 2016). Furthermore, individuals may have sufficient plasticity to adjust during their life-time to withstand different conditions (acclimation). Therefore, future ranges may depend not only on the original species thermal niche (Dell et al., 2011), but also on the possibility of niche shifts by individual adjustments (local acclimation; Peck et al., 2014;Rohr et al., 2018) or, longer-term, adaptive evolution in populations (adaptation) (Conover et al., 2006;Sanford and Kelly, 2011;King et al., 2018) as opposed to niche conservatism across populations (Peterson et al., 1999;Khaliq et al., 2015). Thus, realized thermal niches (based on occurrence data) might vary due to individual modifications through acclimation and are typically smaller than a species fundamental niche based on physiological tolerances, which might be conserved or vary across populations through adaptation. In the case of biological invasions, there is evidence that rapid adaptation (niche shift) can be an important contributor for successful establishment of exotic species (Lee and Bell, 1999;Prentis et al., 2008;Davidson et al., 2011;Parravicini et al., 2015;Wesselmann et al., 2020b), although niche conservatism has also been documented in exotic species (Chefaoui and Varela-Álvarez, 2018;Liu et al., 2020). Therefore, species distribution models (SDMs) predicting the future expansion of exotic species should integrate not only the fundamental niche (physiological tolerance) into correlative SDMs (realized niche) but also include the variation of the fundamental niche across populations in order to consider the potential for niche shifts following invasions (Guisan et al., 2014;Chefaoui and Varela-Álvarez, 2018).
The small tropical seagrass Halophila stipulacea is native to the Red Sea and Indian Ocean (Forskal, 1775;Den Hartog, 1970) and entered the Mediterranean Sea through the Suez Canal over 100 years ago (Fritsch, 1895;Forti, 1927). Currently, H. stipulacea has spread across the eastern (Lipkin, 1975a) and central Mediterranean Sea basin (Sghaier et al., 2011;Gambi et al., 2018). In 2002, H. stipulacea was introduced in the eastern Caribbean Islands probably by anchors of recreational boats (Ruiz and Ballantine, 2004) and recently it reached the northern coast of South America (Vera et al., 2014). This tropical seagrass is able to expand through regions with cooler thermal regimes than those in the native range, which might confer an advantage to invade new environments: in the Mediterranean, H. stipulacea is exposed to much cooler temperatures (minimum SSTs; 14.96 • C) than in the Red Sea (minimum SST: 21.82 • C) (Wesselmann et al., 2020b). As a result, the fundamental thermal niche of H. stipulacea is not conserved within the species geographical range, with significantly lower thermal limits for growth and for metabolism (experimentally determined) for exotic Mediterranean populations than for a native Red Sea population (Wesselmann et al., 2020b). Plants from native and exotic populations were exposed to 12 temperature treatments from 8 to 40 • C and specimens from exotic populations started to show rhizome growth at 14 • C while specimens from a native population at 17 • C (Wesselmann et al., 2020a). This change in the thermal limits occurred in response to the lower Mediterranean SST, as Red Sea populations are not under selective pressure to adapt to lower temperatures, suggesting a rapid adaptation to Mediterranean waters (Wesselmann et al., 2020b). This adaptation to cooler thermal regimes and the trajectories of increase of annual minimum and maximum sea surface temperature (SST) projected for the end of this century (Intergovernmental Panel on Climate Change [IPCC], 2014) suggest that H. stipulacea may be able to spread further and improve its performance in the future due to its high optimum temperature and upper thermal threshold (Wesselmann et al., 2020b).
Here we project the future global distribution of H. stipulacea under projected climatic scenarios using two hybrid models by integrating the experimentally determined thermal limits (fundamental niche), and their variation across native and exotic populations, into the realized species distributions (correlative SDM). Hybrid model 1 was constructed by filtering the correlative SDMs with the physiological thermal limits at the native range (17-36 • C) and Hybrid model 2 by filtering the correlative SDMs with the physiological thermal limits at the exotic range (14-36 • C). These hybrid models, considering niche conservatism (hybrid model 1) or local adaptation (hybrid model 2), were compared to assess differences in H. stipulacea distribution and corroborate or reject the thermal niche shift postulated by Wesselmann et al. (2020b). We hypothesize that the inclusion of the fundamental niche (thermal limits) considering both native and exotic physiological tolerances into the realized species' distributions (correlative SDM) may provide a more robust projection of the spread of this marine exotic species compared to models based on native thermal niches alone.

MATERIALS AND METHODS
Occurrence records of H. stipulacea from its native (Indian Ocean, Red Sea and Arabian Gulf) and exotic range (Mediterranean Sea and the Caribbean islands and northern coast of South America) were extracted from the scientific literature and from online data repositories: the Global Biodiversity Information Facility (GBIF) and the Ocean Biogeographic Information System (OBIF) (Figure 1 and Supplementary Appendix 1). We used the 30 arc-seconds resolution GEBCO's gridded bathymetric data set 1 to delimit the pixels within a 50 m depth, the reported H. stipulacea depth limit (Sharon et al., 2011), as the study coastal area. For georeferencing of H. stipulacea occurrence data and climatic variables, a resolution of 5 arc minutes (∼ 9.2 km) was used, obtaining 553 occurrence records of H. stipulacea entire global distribution, with observations in the exotic range not restricted to harbor, ports and marinas. Only 1 https://www.gebco.net/ FIGURE 1 | Minimum SST and records of Halophila stipulacea included in the current analysis showing its distribution in its native range in the Indian Ocean and Red Sea (blue points) and exotic Mediterranean and Caribbean ranges (pink points). Minimum SST was extracted from bio-oracle.
one record per cell was used to reduce spatial autocorrelation (n = 278 cells).
From those ocean variables available both for current and future conditions, we selected sea surface temperature as the most relevant for limiting the species geographic distribution, according to previous studies (Hemminga and Duarte, 2000;Wesselmann et al., 2020b). Salinity was not included as H. stipulacea is an euryhaline species (Den Hartog, 1970;Hemminga and Duarte, 2000;Oscar et al., 2018), implying that salinity regimes are not expected to limit its spread. Parameters such as fetch or current velocity were also not included in the model because H. stipulacea's distribution is not restricted solely to sheltered areas (Supplementary Appendix 1). Maximum and minimum SST were obtained from Bio-ORACLE v2.0 dataset (Tyberghein et al., 2012;Assis et al., 2018) and tested by Pearson's correlation coefficient < | 0.75| in our study area.
Projected SST values were ensembled from five ocean general circulation models (see Table 1 in  pertaining to the Coupled Model Intercomparison Project Phase 5 (CMIP5 2 ). We used two representative concentration pathways (RCP), the RCP 2.6 and the RCP 8.5, providing the lowest and highest greenhouse gas concentration scenarios, respectively (Moss et al., 2008;Meinshausen et al., 2011). The monthly averaged data by 2050 (from 2030 to 2050) and 2100 (from 2080 to 2100) were computed for each scenario.

Current and Future Species Distribution Models
We implemented an ensemble approach to obtain the realized current and projected distribution of the species by 2050 and 2100 under the RCP 2.6 and the RCP 8.5 scenarios. We used the "biomod2" package (Thuiller et al., 2020) to perform six presence-absence algorithms: flexible discriminant analysis (FDA), generalized additive model (GAM), generalized boosting model (GBM), generalized linear model (GLM), multiple adaptive regression splines (MARS), and randomForest (RF). We extracted at random the same number of pseudo-absences as presences for each of the 3 sets of pseudo-absences used. Data were split into a calibration (70%) and a validation (30%) set in each of the 5 iterations performed for each model and set of pseudo-absences. Thus, 90 models were computed (3 pseudo-absence sets × 6 modeling techniques × 5 iterations). We evaluated model performance using the true skill statistic (TSS; Allouche et al., 2006), the area under the receiver operating characteristic (ROC) curve (AUC), and ROC-derived sensitivity (presences correctly projected) and specificity (absences correctly projected) measures (Fielding and Bell, 1997), with the threshold which optimized ROC and TSS scores (Thuiller et al., 2020). To produce the ensemble, we just used the models which obtained TSS > 0.6. The consensus projection under current conditions was an ensemble computed through the average of binary projection (committee averaging ensemble). Ensembles were projected to the RCP 2.6 and 8.5 scenarios to obtain the projection for the future. To assess the uncertainty in future projections related to the different values for each variable between the training range and the novel scenarios we computed a clamping mask.
We estimated the importance of each variable following a procedure similar to "randomForest" (Liaw and Wiener, 2002). The correlation between the full model and a model rearranged without one variable using three iterations was calculated to obtain an importance value from 0 to 1 (highest importance) (Thuiller et al., 2020). Analyses were performed in R (R Core Team, 2020) using a computational cluster facility. We applied the true skill statistic (TSS) related threshold to transform projection into binary maps.

Hybrid Models
Halophila stipulacea is able to survive in regions with SST from 8 to 36 • C but the species cannot grow properly and develop viable populations at the lower range of this thermal interval for survival (Wesselmann et al., 2020a). Thus, we used the observed thermal limit of rhizome growth for H. stipulacea reported in the dataset of Wesselmann et al. (2020a) for the Red Sea (SST = 17-36 • C) and the Mediterranean Sea (SST = 14-36 • C) to produce the hybrid models, where the lower thermal limit for rhizome growth in Mediterranean populations compared to those in the Red Sea provides evidence of adaptation when growing in cooler, exotic ranges.
We applied these thresholds to reclassify the ensembles of minimum and maximum SST pertaining to each RCP scenario by 2050 and 2100 using the "raster" package (Hijmans, 2020) in R. Finally, we filtered the projection obtained by the correlative SDM (realized niche) with the SST related thresholds of growth of native and exotic populations (fundamental niche) to obtain two hybrid models: i) Hybrid model 1 was obtained by filtering the correlative SDM with the thermal limits of growth reported in the native range (SST = 17-36 • C), considering niche conservatism and ii) Hybrid model 2 by filtering the correlative SDM with the limits of growth reported for the exotic Mediterranean populations (SST = 14-36 • C), allowing, therefore, for local adaptation based on experimental evidence on performance for Mediterranean populations (Wesselmann et al., 2020b).
For each hybrid model, we estimated the areas of agreement (area identified as suitable by the correlative model and the thermal limits, i.e., where the realized and fundamental distribution match), overestimation (area identified as suitable by the correlative model but with SST outside the thermal limits for growth of the species) and underestimation (area not identified as suitable by the correlative model but with SST inside the thermal limits for growth of H. stipulacea) of the correlative SDMs in relation to the thermal limits of growth as a percentage of cells.

RESULTS
The highest validation scores of the models were produced by Random Forest (RF) with TSS (± SD) = 0.567 ± 0.074 and the lowest by Generalized Linear Model (GLM) with TSS = 0.39 ± 0.068 (Supplementary Table 1). The mean evaluation scores for the ensemble were: AUC = 0.957, sensitivity = 87.008, and specificity = 94.712 and TSS = 0.817. Minimum sea surface temperature was the variable with the highest relative importance (mean scores = 0.868) for the ensemble model, compared to the maximum sea surface temperature (0.66). This adds to previous evidence supporting minimum SST a limiting factor of the current distribution of H. stipulacea in the Mediterranean Sea (Wesselmann et al., 2020b).

Current and Future Correlative Species Distribution Models
The ensemble model for the present time approximated the actual distribution of the species, matching a high proportion of the reported species occurrences in the native range (78%) and in the exotic Caribbean (90%) and Mediterranean (82%) ranges (Supplementary Figure 1). Min SST extracted from the thermal occurrence envelopes used in the correlative model was lower for the exotic Mediterranean than for the exotic Caribbean and the native range (Figure 2 and Supplementary Table 2). The model found suitable the Indian Ocean, Red Sea, Arabian Gulf, eastern Caribbean Islands and eastern Mediterranean Sea, as well as most of the western Mediterranean, coast of Mexico, Florida, South Carolina and the western Caribbean Islands (Supplementary Figure 1), where the species has not been reported so far (Figure 1). The correlative model projected future contractions of H. stipulacea suitable area in the north of the Red Sea, Arabian Gulf, Caribbean Islands and eastern Mediterranean by 2050 (RCP 2.6 and 8.5) and 2100 (RCP 2.6) and the complete disappearance from the eastern Mediterranean under the most severe scenario by 2100 (under RCP 8.5; Supplementary  Figure 1). Habitat increases were projected for the south of the Red Sea, Seychelles, western Mediterranean and Atlantic coast of Morocco by 2050 and 2100 (RCP 2.6). The ensemble also identified the western coast of Africa as presently suitable for H. stipulacea (Supplementary Figure 1). Hence, the currently realized niche is much smaller than the potential suitable niche. The results suggest loss from the Red Sea, which would become too warm to support the species. However, clamping masks revealed uncertainty in some areas of the Red Sea due to the different values of one variable between the training range and the worst greenhouse gas emission scenario (RCP 8.5 by 2100; see Supplementary Figure 2). Thus, future projections for this region and scenario should be taken with caution.

Spatial Distribution of Thermal Limits for Growth
The lower thermal limits would restrict H. stipulacea distribution, whereas upper thermal thresholds would not constrain the predicted distribution of the species in the future, as all coastal cells would present maximum SST below 36 • C, even under the most severe scenario in 2100. The lower thermal limits for growth of H. stipulacea in its native range implied, under present conditions, suitable temperatures for growth in the western Pacific from 15 • N to 23 • S, in the western Atlantic from 35 • N to 28 • S and in the eastern Atlantic through the coast of Africa (Supplementary Figure 1). The lower thermal limits of H. stipulacea obtained from Mediterranean populations also identified as suitable the northern coast of Argentina and North Carolina (Supplementary Figure 1). This suitable area for growth would increase in the western Atlantic by 2100, from North Carolina to Virginia (United States) under RCP 2.6 and to Pennsylvania (United States) under RCP 8.5 by 2100 (Supplementary Figure 1).
In the Mediterranean Sea, according to the thermal limits for growth of the species in its native range (

Hybrid Model
Overall, the output from our hybrid models revealed that the correlative model underestimated H. stipulacea potential areal extent in many regions, as maximum seawater temperatures will not exceed the species thermal limits (36 • C) in any future scenario. Therefore, H. stipulacea would not lose habitat and current seagrass meadows would persist in the future, maintaining 99% of the cells occupied at present (Figure 3). The hybrid model based on native Red Sea thermal limits and on exotic Mediterranean thermal limits (Hybrid model 1 and 2, respectively, hereafter) provided identical results for the Indian Ocean, relatively similar distributions for the Caribbean Sea and projections diverged for the Mediterranean Sea (Figure 3). The area of agreement and underestimation of Hybrid model 1 and 2 in the present distribution identified most of the reported species occurrences within the native (87%) and Caribbean exotic ranges (97%), exceeding the outcome of the correlative SDMs by 10% (Figures 1, 4). In contrast, in the Mediterranean Sea, the hybrid model 1 (agreement and underestimation area) recognized solely 13% of the reported occurrences of H. stipulacea, while the hybrid model 2 identified 90% of the occurrences (Figures 1, 5).
In the native range, both hybrid models projected an increase in the area of agreement from the present (13%) to 2050 (16-17%) and to 2100 under RCP 8.5 (22%) and identified an extensive area underestimated by the correlative model, accounting for 76-83% of the cells along the coast of the native range at present and future scenarios (Figure 3). Considering the area of agreement and underestimation, in all future warming scenarios, H. stipulacea could gain 1-2% of its habitat in the Arabian Gulf (Figures 3, 4 and Supplementary Figure 3).
In the Caribbean, the area of agreement between the physiological thermal limits and the correlative models in the hybrid model 1 increased from 7% in the present to 15% by 2100 under RCP 8.5 scenario, while the agreement area in the hybrid model 2 was wider (from 9 to 18%) because it included the coast of Argentina and the coast of North Carolina as suitable to support H. stipulacea's growth (Figures 3, 4 and  Supplementary Figure 3). The correlative model underestimated a huge potentially suitable area, corresponding to 65-74% of all total cells of the western Atlantic coast at the present and under future RCP scenarios, according to hybrid model 1 and 2 (Figures 3, 5 and Supplementary Figure 3). Taking into account the agreement and underestimated area, H. stipulacea could spread southwards and westwards, reaching the Dominican Republic, Cuba, the Gulf of Mexico, Venezuela and Brazil, gaining 7% of its habitat in 2050 and 11% in 2100 under the most severe scenario (RCP 8.5) according to hybrid model 1, whereas according to hybrid model 2, it could increase from 3 to 8% (Figures 3, 4 and Supplementary Figure 3).
In the Mediterranean Sea, the hybrid model 1 showed a reduced suitable area for H. stipulacea in comparison with the hybrid model 2. The area of agreement in Hybrid model 1 was negligible at present and under all future warming scenarios (0-3% of all the cells) (Figures 3, 5) and for most of the western Mediterranean Sea it was overestimated by the correlative model and identified as unsuitable for H. stipulacea's growth, as minimum temperatures are below 17 • C (Figures 3, 5). Even adding the area of agreement and underestimation, the suitable habitat would still be restricted to 4% of all the cells of the Mediterranean at present, 10-12% by 2050 and 2100 (RCP 2.6), and 33% by 2100 (RCP 8.5, Figure 3, 5). Conversely, the agreement area of hybrid model 2 at present, 2050 and 2100 (RCP 2.6) occupied most of the eastern Mediterranean basin and the south-western basin along the coast of Africa (25-35% of all the Mediterranean coastal cells) and disagreed with the correlative model prediction of the presence of the species at the southeastern Iberian Peninsula, southern Adriatic and northern Aegean Sea, as temperature at these regions would remain unsuitable for H. stipulacea growth (< 14 • C). By 2100 under RCP 8.5, H. stipulacea could spread into the western Mediterranean Sea, extending from Valencia to Andalusia and the Balearic Islands, and into the eastern Atlantic Ocean through the Gibraltar strait, although the area of agreement of hybrid model 2 would decrease to 5% of all the cells in the Mediterranean Sea (Figures 3, 5). This reduction is due to an underestimation of the correlative model, as 15% at present, 22-24% in 2050 (RCP 2.6 and 8.5) and 2100 (RCP 2.6) and up to 72% by 2100 (RCP8.5) of all the cells in the Mediterranean Sea, mostly concentrated in the eastern basin, would be suitable for H. stipulacea according to the thermal limits of growth of the species (14-36 • C, Figure 5). Therefore, adding together the area of agreement and underestimation of hybrid model 2, H. stipulacea could actually occupy 50% of all the cells of the Mediterranean coast at present and increase its habitat by 51% by 2100 under RCP 8.5, occupying 76% of all the cells of the Mediterranean spreading through the western basin (Figures 3, 5).

DISCUSSION
In this study, we show that the global future spread of Halophila stipulacea in contrasting climate change scenarios is better predicted when accounting for local adaptation as compared to niche conservatism. This novel outcome was possible by  combining correlative SDM with experimentally resolved estimates of the thermal limits for growth of Halophila stipulacea, obtained for the exotic Mediterranean range and for the native range in the Red Sea. The hybrid models derived from correlative models and thermal niches of populations at the native range (hybrid model 1: niche conservatism) and the Mediterranean range (local adaptation) were similar for the native range (Indian Ocean and Red Sea) and exotic Caribbean range, whereas in the Mediterranean Sea the hybrid model based on the thermal limits of the native range failed to predict most of the current occurrences of the species, which were however, correctly captured by the model based on the thermal limits of the exotic populations in the Mediterranean range that incorporated local adaptation. This study demonstrates the importance of including knowledge on the physiological thermal limits of the species and accounting for their variation across populations, and adaptation and acclimation capacity when predicting future distributions, particularly for biological invasions at regions where thermal conditions extend beyond the limits of the native range.
Hybrid models projected that H. stipulacea will maintain all its current habitat in the future and increase it slightly in the exotic Caribbean (7% by 2050 and 11% by 2100) and up to 50% in the Mediterranean ranges, expanding through the western Mediterranean basin (except the Gulf of Lion) under the most severe scenario (RCP 8.5) by 2100, as previously hypothesized for the species (Wesselmann et al., 2020b). Similar westward longitudinal expansions in response to climate-driven temperature increases have been projected for other Lessepsian migrants, such as fish (Pais et al., 2007;Azzurro et al., 2013Azzurro et al., , 2017Coro et al., 2018;Amen and Azzurro, 2020), foraminifera (Langer et al., 2012) and bivalves (Sarà et al., 2018).
The current distribution of H. stipulacea in the Mediterranean Sea is related to the capacity of the Mediterranean populations to grow at 14 • C (lower than the native populations; Wesselmann et al., 2020b), allowing its westward expansion. However, the low winter temperatures still prevent occurrence of H. stipulacea in most of the western Mediterranean basin at present. Local adaptation has been described as a mechanism likely employed by exotic species that establish successfully in new environments (Prentis et al., 2008). However, this process has been ignored by niche prediction models, largely because of a gap of experimental data on local adaptation of the thermal niche of exotic species (Valladares et al., 2014).
Our results also show that the expansion of H. stipulacea beyond its native range is not only limited by its capacity to tolerate certain temperatures but is also dispersal limited. Our distribution models of H. stipulacea indicate the suitability of multiple regions beyond the range of current occurrences (Figures 1, 3), such as the coast of North Carolina, the Caribbean coast of Honduras and Nicaragua and the Algerian coast in the western Mediterranean. This suggests that the occupation of potential areas with suitable temperatures is limited by long distance dispersal, which could occur through rafting of rhizome fragments (Steiner and Willette, 2015;Smulders et al., 2017;Willette et al., 2020) or seeds (Chiquillo et al., 2018;Nguyen et al., 2018), despite possible constraints within short scales by locally acting factors (e.g., seawater turbidity and exposure). Long-range seed transport seems unlikely as the genera Halophila produces dormant seeds, neutral or negatively buoyant, and are therefore transported largely through bedload transport over short distances (up to 100 of meters; Kendrick et al., 2012), although longer ranges may be reached during storms or if ingested by highly mobile grazers (Heck et al., 2008;Tol et al., 2017). Successful long-term dispersal through rhizome fragments appears viable, as small (2 cm) fragments are able to grow floating in the seawater in the laboratory  and a field experiment demonstrated settlement of H. stipulacea fragments within 10 days of detachment (Smulders et al., 2017). Fragments could reach considerable distances over 10 days, either drifting or inadvertently transported by vessels or other vectors. Even a single fragment could create a new clonal population and disperse even further through a stepping-stone model, as hypothesized for other seagrasses (Arnaud-Haond et al., 2012). The prevailing Caribbean Current could facilitate the colonization of the western Caribbean, which is in agreement with a prediction made by Ruiz et al. (2017). However, the further colonization of the coast of Venezuela and the coast of Brazil through currents is less likely due the westward direction of the Caribbean Current (Richardson, 2005). Similarly, considering the water circulation in the Mediterranean Sea, H. stipulacea could colonize first the entire western Italian coast and Corsica than the Algerian coast because the oceanographic current flows counterclockwise in the western Mediterranean (Millot and Taupier-Letage, 2005). Despite this, human mediated transport would be required to fill the potential niche of H. stipulacea populations, as shown by its arrival to Martinique and the Lesser Antilles through recreational boats (Ruiz and Ballantine, 2004;Willette and Ambrose, 2009) or to other H. stipulacea populations in the Mediterranean dispersed either as fragments in fishing gear or attached to the hull of ships or boat anchors (Lipkin, 1975b;Gambi et al., 2018).
Although the current presence of H. stipulacea in the Mediterranean coast is limited by the lower sea temperatures (Wesselmann et al., 2020b), the present study supports the hypothesis that H. stipulacea will continue to expand into the western basin propelled by increasing SST (Wesselmann et al., 2020b). However, our results differ from previous studies predicting a future expansion of H. stipulacea along the northwestern Mediterranean coast already by 2050 (Beca-Carretero et al., 2020) or by 2100 (Gamliel et al., 2020). The capacity of H. stipulacea to both cope with very warm waters and to grow at temperatures several degrees below the minimum temperatures in its native range are key to successfully continue to spread. Exotic marine species often have a greater thermotolerance than related native species (Sorte et al., 2010;Zerebecki and Sorte, 2011;Magozzi and Calosi, 2015), and can profit from warming by expanding in their recipient range. Indeed, the upper thermal tolerance of H. stipulacea is higher than native Mediterranean seagrasses, especially for Posidonia oceanica (Savva et al., 2018). Future projections indicate temperature increases that would exceed the physiological upper thermal tolerance limit of P. oceanica and lead to a decline of its habitat (Jordà et al., 2012;Marbà et al., 2015). Recently, it has been projected that, under the most severe climate change scenario, C. nodosa could lose 46% of its habitat and P. oceanica's habitat could be completely lost by 2100 . If these projections are met and the native plants do not adapt, the Mediterranean will lose key foundation species, which may be replaced by H. stipulacea. The future expansion of H. stipulacea and the potential shift in the dominant flora in the Mediterranean (from P. oceanica to H. stipulacea and the native warm tolerant species C. nodosa) may have profound effects on seagrass ecosystem functioning. The ecological roles of P. oceanica and H. sipulacea are different as they have different evolutionary histories and H. stipulacea tends to colonize muddy estuarine waters and unstable environments (Wesselmann et al., 2021). The limited information available on the provision of ecosystem functions of H. stipulacea in the Mediterranean Sea indicates a similar carbon storage capacity of H. stipulacea and native seagrasses in the eastern Mediterranean (Wesselmann et al., 2021), although carbon stocks in that region ranks among the lowest across seagrass meadows (Fourqurean et al., 2012;Duarte et al., 2013;Mazarrasa et al., 2017). Besides, the dense and tall canopy of P. oceanica attenuates the waves and its meadows play a key role on coastal protection (Sánchez-González et al., 2011), whereas the provision of this ecosystem service by H. stipulacea, with a canopy few centimeters high (Den Hartog, 1970), is low (Christianen et al., 2013). Moreover, the three-dimensional structure of P. oceanica canopy provides food and more shelter for diverse faunal assemblages (Francour, 1997;Boudouresque et al., 2006) than the short one of H. stipulacea. Faunal assemblages in H. stipulacea meadows are mainly characterized by species from the phyla Mollusca, Amphipoda, Decapoda, Polychaeta, Crustacea (Cancemi et al., 1994;Acunto et al., 1997) and a fish community dominated by sparids, labrids, and gobiids in which planktivorous species (common in other Mediterranean seagrasses) are absent (Di Martino et al., 2007). Di Genio et al. (2021) found higher species richness of benthic fishes in H. stipulacea meadows than in the nearby dead matte of P. oceanica, however, it remains unclear if the fish abundance and species richness is similar in H. stipulacea and P. oceanica meadows. Finally, H. stipulacea may also drive populations shifts in species that utilize P. oceanica as a component of their diet (Tomas et al., 2005), although it has already been documented that the herbivorous fish Sarpa salpa, typically feeding on P. oceanica, grazes year-round on H. stipulacea leaves (Di Genio et al., 2021) and signs of grazing were recorded in more than half of all H. stipulacea leaves (Gambi et al., 2018).
On the contrary, in the exotic range in the Atlantic Ocean, H. stipulacea does not show higher thermal tolerance than native seagrasses. While in the Mediterranean H. stipulacea coexists with temperate seagrasses (although C. nodosa has a warm affinity; Savva et al., 2018), in the native range and in the Caribbean Sea it coexists with native, warm tolerant species (Den Hartog, 1970;Willette and Ambrose, 2012). As tropical seagrasses have similar thermal tolerances (McMillan, 1984;Koch et al., 2007), sustained elevated temperatures in the Caribbean Sea are unlikely to promote the replacement of the native seagrasses by H. stipulacea. In the Caribbean Sea, H. stipulacea mostly colonized bare "sand halos, " where it forms monospecific meadows, although it has also been observed in mixed seagrass beds with Thalassia testudinum, Syringodium filiforme, Halodule wrightii, and H. decipiens (Willette et al., 2014). Therefore, in this region the expansion of H. stipulacea might be driven by its fast growth rates (Wesselmann et al., 2020b), low light requirements (Sharon et al., 2011), capacity to resist siltation (Terrados et al., 1998) and fast recovery after perturbations (e.g., turtle grazing, storms, anchoring, or bioturbation: Steiner and Willette, 2015;Smulders et al., 2017;Christianen et al., 2019).
In general, the distribution of H. stipulacea based on thermal limits delimited a wider potential range for the species than the realized native and exotic niche. Our hybrid models, combining correlative parametrization with experimentally derived evidence of thermal performance in native and exotic ranges, clearly demonstrate that overlapping thermal limits with SDM is useful to provide an improved future projection on responses of exotic species to climate warming, as it is able to address regions where statistical approaches mismatch physiological information. Thermal limits have been combined before with SDM to improve projections (e.g., Woodin et al., 2013;Martínez et al., 2015;Talluto et al., 2016;Gamliel et al., 2020). However, the thermal tolerance of species over a wide thermal range and its potential variation due to local adaptation are usually ignored in these models (Valladares et al., 2014). Gamliel et al. (2020) combined SDMs with data on H. stipulacea's physiological thermal limits, obtained after exposing plants to a narrower temperature range (from 10 to 30 • C; Georgiou et al., 2016) than those used here (from 8 to 40 • C; Wesselmann et al., 2020a), to project H. stipulacea expansion in the northern Mediterranean coast. Our study takes the additional step of combining SDM with knowledge on thermal limits, considering niche conservatism and local adaptation, to project future trajectories of marine invasions. However, the thermal limits for rhizome growth used in the hybrid model do not consider the phenology of the species and that local populations could potentially adjust their rhizome growth to match the optimal conditions for growth, avoiding the lower thermal limit. Moreover, the integration of other relevant parameters, such as reproductive phenology (Chefaoui et al., 2019), oceanographic currents (Assis et al., 2015), or biotic interactions (Record et al., 2018) would certainly improve projections.

CONCLUSION
In conclusion, including the variation of thermal limits across native and exotic populations in SDMs improves the prediction of realized, and hence likely projected, distributions of exotic species. Our models support two important insights, (1) warming threatens the expansion of H. stipulacea habitat suitability through the western Mediterranean under the most severe scenario (RCP 8.5) by 2100, increasing its distribution by 50% in the Mediterranean; (2) the potential range of H. stipulacea is vastly larger than that presently occupied, so risks of contemporary expansions in subtropical and tropical waters are high. In the future, H. stipulacea may replace, but not displace, dominant temperate seagrasses in the Mediterranean Sea, as the thermal limits of the native species are exceeded with future climate warming, resulting in important changes in the ecosystem functioning of Mediterranean seagrass meadows .

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
CD, NM, ES, and MW conceived the study design. RC developed the models. MW produced the figures of the models and analyzed the data. CD, NM, ES, RC, and MW interpreted the results. MW wrote the initial draft and all coauthors contributed to editing the manuscript and gave final approval for publication.

FUNDING
This study was funded by the Spanish Ministry of Science, Innovation and Universities (SuMaEco RTI2018-095441-B-C21), King Abdullah University for Science and Technology (3834 KAUST-CSIC Research Collaboration), the Foundation for Science and Technology (UIDB/04326/2020), and the Operational Programmes CRESC Algarve 2020 and COMPETE 2020 (EMBRC.PT ALG-01-0145-FEDER-022121 and BIODATA.PT ALG-01-0145-FEDER-022231). MW was supported by a Ph.D. contract (BES-2016-078241) of the Spanish Ministry of Science, Innovation and Universities. RC obtained additional support from the European Union's Horizon.

ACKNOWLEDGMENTS
We thank Cymon Cox for his help with the CCMAR highperformance computing infrastructure (CETA).