Invasion and Extirpation Potential of Native and Invasive Spartina Species Under Climate Change

Coastal areas host some of the planet’s most productive ecosystems, providing life-sustaining ecological services and several benefits to humankind, while also being some of the most threatened areas (e.g., by globalization, climate change, and biological invasion). Salt marshes are coastal habitats with a key role in food and shelter provisioning, sediment deposition, nutrient cycling and carbon storage. Spartina spp. is a genus of grass halophytes which occurs in salt marshes worldwide, and includes species with different invasive potential. We evaluated the effect of climate change in the distribution and invasion potential of five Spartina species (S. anglica, S. alterniflora, S. densiflora, S. patens, and S. maritima) at a global scale. Species distribution models (SDMs) were applied on species occurrence data and atmospheric environmental predictors (WorldClim 2.1) to project potential changes in habitat suitability and associated changes in distribution and species co-occurrence until the end of the century, across four Shared Socioeconomic Pathway scenarios (i.e., SSP1-2.6 to SSP5-8.5). Projections showed a global trend for increasing species co-occurrence, with a general range expansion potentiated by increasing pathway severity. This study suggests that Spartina species can potentially benefit from climate change, predicting poleward expansions in the Northern Hemisphere for most species, with results pointing at increased conflict and invasion potential in Northern Europe and East Asian shorelines, already under strong invasive pressure. S. anglica is projected to remain a successful invader, with more severe scenarios likely favoring greater expansions. S. alterniflora exhibits very low expansion comparatively, despite exhibiting the same northward distribution shift. SSP1-2.6 produced the smallest change to species co-occurrence, suggesting a smaller potential for invasion-related conflicts, although still registering a potential net expansion for the Genus. Despite their limitations, SDMs can help establish general trends in climate change ecology and inform policymakers and environmental agents to ensure the correct management of these habitats and, ultimately, ecosystems.


INTRODUCTION
Coastal zones present some of the planet's most productive ecosystems, even though they occupy a relatively small percentage of total land and ocean extension (Reid, 2005;Himes-Cornell et al., 2018). These land-ocean transitional habitats harbor vegetated assemblages that provide a vast array of lifesustaining ecological services and several benefits to humankind (Barbier, 2013;Himes-Cornell et al., 2018), and include salt marshes, mangroves, and seagrass beds, which are known to play a key role in large scale biogeochemical processes such as carbon sequestration and nutrient cycling (Duarte et al., 2005(Duarte et al., , 2021Mueller et al., 2019). Considering the proximity of these ecosystems to human-related activities, it is not surprising that coastal areas include some of the world's most threatened ecosystems, with studies showing that nearly half of salt marshes, 35% of mangroves and 29% of seagrass beds have been lost in the past 50 years alone (Himes-Cornell et al., 2018). Salt marshes are distributed worldwide, occurring from polar to tropical regions (Mcowen et al., 2017), and establish in lowenergy intertidal zones located on the fringes of the inner areas of bays and estuaries . This habitat type typically consists of assemblages between halophytic (i.e., salttolerant) herbs, grasses, and low shrubs, which are adapted to regular or occasional tidal immersion, as well as a great diversity of land and marine wildlife (Mcowen et al., 2017). It also plays a key role providing a wide array of ecosystem services such as habitat provision to migratory birds and vast plant and animal assemblages (Greenberg et al., 2014); sequestration of carbon and other pollutants (Duarte et al., 2010(Duarte et al., , 2021; storage of sediments and nutrients (Duarte et al., , 2014; water purification, storm surge, and flood attenuation (Zedler and Kercher, 2005;Donatelli et al., 2018;Leonardi et al., 2018).
The proximity of these ecosystems to urban settlements means that salt marshes are highly vulnerable to increased pressure from anthropogenic sources (Mcowen et al., 2017;Tonti et al., 2018;Curado et al., 2020). Due to global warming, high rates of environmental degradation, and rising sea levels, together with encroachment by invasive species have increasingly being reported (Duarte et al., 2015;Himes-Cornell et al., 2018;Curado et al., 2020). Indeed, salt marshes and other coastal wetlands are among some of the ecosystem types most affected by biological invasions (Gedan et al., 2009;Vilà et al., 2011). The spread of invasive species can have broad-scale impacts in these areas, affecting community dynamics and physically changing the marsh structure, potentially reducing biological diversity, and permanently altering their ecological functions and services (Gedan et al., 2009).
Cordgrasses (of the Genus Spartina, Poaceae), are a group of halophyte species widespread in salt marshes around the world (Castillo et al., 2017;Bortolus et al., 2019). These grasses colonize intertidal mudflats and are known to be important ecosystem engineers due to their ability to stabilize sediments and contribute to salt marsh formation (van Hulzen et al., 2007;Vu et al., 2017). Indeed, this group is of great importance, since these species facilitate ecological succession, and provide marshes with the spatial structure needed to promote high biodiversity levels (Castillo et al., 2017). However, cordgrasses are also prolific invasive halophytes (Strong and Ayres, 2013;Castillo et al., 2017). Due to their physiological characteristics [e.g., root systems adapted to soil inundation and low oxygen levels (Castillo et al., 2017)], species of the Genus Spartina can easily dominate frequently flooded tidal flats (Pennings et al., 2005). Additionally, cordgrasses propagate through both sexual and clonal reproduction. While clonal growth strategy partially determines the invasive potential of this group of halophyte plants (Castillo et al., 2016), seeds, and rhizome propagules differing in size, resource requirements, and genetic identity result in different invasive potential across species (Castillo et al., 2017). Together with biomass accumulation, both aspects are responsible for invasive Spartina species overgrowing native salt marshes, leading to significant decreases in plant and fauna abundance and diversity, as well as carrying the potential to hybridize with native Spartina species (Strong and Ayres, 2013). Moreover, the presence of these invasive Spartina species is also recognized to have severe impacts on the salt marsh ecosystem services, such as its important role in the natural remediation of transitional ecosystems (Human et al., 2020).
Some areas are known to be heavily invaded by cordgrass [e.g., the Chinese coast (Strong and Ayres, 2013)], with records of these invasions going back several decades, or even centuries (Castillo et al., 2018). Some of the most studied invasive cordgrass species include the common cordgrass, Spartina anglica (Hubb); the smooth cordgrass, Spartina alterniflora (Loisel); the austral cordgrass, Spartina densiflora Brongn; and the salt meadow cordgrass, Spartina patens (Aiton) Muhl. The common cordgrass, S. anglica, is a highly aggressive non-native species (Strong and Ayres, 2013). This species emerged through allopolyploid speciation (Guenegou et al., 1988;Ainouche et al., 2004), originating from the known hybrid S. × townsendii -the F1 hybrid of Spartina maritima and S. alterniflora (Strong and Ayres, 2013) and is endemic to the British Islands. Despite having been shown to have lower growth performance and competitive ability when compared to other invasive cordgrasses (such as S. alterniflora) in some areas, S. anglica has spread from the coasts of Europe to North and South American coastlines, to China, and even to the Australian and New Zealand shore (Bortolus et al., 2019),being responsible for some of the largest continentscale invasion events in Asia and in the Americas, with major ecological and socioeconomic impacts (Bortolus et al., 2019).
The smooth cordgrass, S. alterniflora, is present in temperate and subtropical salt marsh habitats (Bortolus et al., 2019), being native to the Atlantic and Gulf coasts of North and South America (Kirwan et al., 2009) and dominating intertidal marshes (Zheng et al., 2018). Like S. anglica, S. alterniflora is also an invasive cordgrass, both inadvertently and deliberately introduced in coastal wetlands worldwide (Zhang et al., 2017;Bortolus et al., 2019). Over the past centuries, reports of invasion of European coastlines, mainly in France and Southern England (Bertness, 1991;Baumel et al., 2003), in the Mediterranean (Hessini et al., 2009), China, andNew Zealand (Feng et al., 2017;Zheng et al., 2018) have been increasing; and S. alterniflora is now very well established in most of these areas due to its high adaptability and propagation ability (Zuo et al., 2012;Zhang et al., 2017). The saltmeadow cordgrass, S. patens, is a native cordgrass to brackish salt marshes and coastal dunes of the Atlantic and Gulf coasts of North America (Bertness, 1991). This species has now a vast non-indigenous distribution in coastlines of Europe, along the western Mediterranean coast and the Atlantic coast of the Iberian Peninsula (Baumel et al., 2016;Duarte et al., 2018). Finally, the denseflower cordgrass, S. densiflora, is native to the South American continent and has spread to Europe, in particular to the region of the southwestern Iberian Peninsula (Nieva et al., 2001;Bortolus, 2006), Morocco, Mexico, and the western coast of the United States (Bortolus, 2006). The species is known to have hybridized with the small cordgrass (S. maritima) in the salt marshes of Southern Iberian Peninsula (Castillo et al., 2018) and is a case study for the potential latency of the invasion process of plant species followed by unprecedented rates of spreading (Duarte et al., 2018;Bortolus et al., 2019). An example of endemic cordgrass species suffering from invasion from other sister species is that of S. maritima. This is a known native of western and southern European shores (from the Netherlands to the south of England and Ireland), and is the only cordgrass species that is native to the salt marshes of continental Europe (Castillo et al., 2018). It also occurs along the Mediterranean Sea coastlines to the Atlantic coast of Morocco, and two separated populations in the Atlantic coasts of Namibia and South Africa (Bortolus et al., 2019). When competing with other cordgrass invasive species, S. maritima is known to have suffered hybridization at least twice: (1) with S. alterniflora, originating S. × townsendii, which later led to S. anglica; and (2) with S. densiflora; which poses an increased threat of replacement in its native habitats.
The potential for climate change together with biological invasions to affect cordgrass species distribution is of particular interest. The saltmeadow cordgrass (S. patens) is rapidly disappearing from wetlands in the coastal marshes of the United States due to current inundation patterns, and accelerated sea-level rise will likely reduce this species' persistence (Watson et al., 2016). However, temperature rise does not seem to significantly impact S. patens or its role as a foundational species (Gedan and Bertness, 2010;Duarte et al., 2016). The hybrid formation has been shown to change in response to changes in atmospheric variables. Specifically, the establishment of reciprocal hybrids between the native S. maritima and the invasive S. densiflora in SW Iberian Peninsula was shown to change with air temperature and rainfall, with a negative effect on the native and endangered species S. maritima (Duarte et al., 2015;Gallego-Tévar et al., 2019). Reciprocal transplant trials also revealed that S. alterniflora might be stressed by future warming, leading to decreased belowground allocation and peat accumulation (Crosby et al., 2017).
An accurate description of a species' ecological and geographical distribution is fundamental for understanding existing patterns of biodiversity, and the processes responsible for shaping them (Wisz et al., 2008). Given the differential effects from future conditions in cordgrass fitness, it is important to predict how cordgrass species might behave under climate change scenarios, in particular regarding changes to their potential distribution. By analyzing predicted distributions, it is possible to recognize where species conflicts may emerge, and identify potential areas with higher potential of invasion. This has become a very important part of conservation planning in the past decades (Elith et al., 2006), with a wide array of modeling techniques being applied in predicting species or community distribution, particularly for plant species (Bagheri et al., 2017). Species distribution models (SDMs) have seen a steep rise in their development and use (Zimmermann et al., 2010). While bounded by a set of assumptions and limitations that should be considered during interpretation (Araújo and Pearson, 2005;Araújo and Guisan, 2006;Heikkinen et al., 2006;Fitzpatrick and Hargrove, 2009), these models are widely employed and offer a good framework for predicting species distributions and their changes, regarding large species assemblages and vast geographical spaces (Elith et al., 2006). This type of models has now been used for several decades to determine suitable habitats for species or communities, namely to evaluate the invasive potential of non-indigenous species and their proliferation (Peterson and Vieglais, 2001;Araújo and Rahbek, 2006;Jiménez-Valverde et al., 2011;García-Roselló et al., 2014) and project the potential effects of climate change on the distribution of species or communities (Araújo and Rahbek, 2006;Peterson and Soberón, 2012;García-Roselló et al., 2014). Notably, SDMs that incorporate future climate change predictions are considered an effective way to address some of the questions behind climate change effects on biodiversity (Sinclair et al., 2010). Given this context, the present study aims to evaluate the potential biogeographical impact of climate change on cordgrasses worldwide. Specifically, this study implements a series of SDM ensembles to estimate how future conditions may affect habitat suitability, species richness and colonization, and extinction susceptibility for five species of the Spartina Genus [S. anglica (Hubb), S. alterniflora (Loisel), S. densiflora Brongn, S. patens (Aiton) Muhl, and S. maritima (Curtis) Fernald] at a global scale. This analysis is based on the 30 year average of projections for atmospheric variables (monthly average maximum temperature; monthly average minimum temperature; and monthly average total precipitation). The analysis encompasses both an historical timeframe (early 2000s) and different future time-spans (i.e., 2021-2040, 2041-2060, 2061-2080, and 2081-2100), considering four of the Shared Socioeconomic Pathways (SSP1-2.6, 2-4.5, 3-7.0, and 5-8.5) of the Coupled Model Intercomparison Project phase 6 (CMIP6) (Tebaldi et al., 2021).

Species Data and Curation
The occurrence data used in this study was compiled from data collected from the Global Biodiversity Information Facility (GBIF), and from literature surveyed from the Web of Knowledge in the case of specific species: S. densiflora, S. maritima, and S. patens in the Iberian Peninsula ; and S. alterniflora in the coastline of China (Cheng et al., 2006;Gao et al., 2007;Liu et al., 2007;Zhang et al., 2010Zhang et al., , 2020Wang et al., 2013Wang et al., , 2015Yang et al., 2013;Jin et al., 2016;Pang et al., 2017;Maebara et al., 2020;Xia et al., 2020;Yue et al., 2021 A specific set of filters was used for all species in each dataset to limit the data retrieved. The applied filters included "Human Observation" and "Preserved specimen" (for the Basis of Record); "Has Coordinate" set as true; "Has Geospatial issue" set as false. The endemic status (i.e., endemic versus non-indigenous species) was also retrieved for each occurrence point. The occurrence dataset, together with the plotted dataset occurrences are present in Supplementary Material (see "Spartina_Rfull.csv" and "plotted_occur.zip, " respectively). The compiled dataset was then curated using R studio. First, the original dataframe was converted into a spatial polygon object, which was then restricted by removing all points outside the extent of the coastline inland-buffer shapefile. This restriction procedure of both species' occurrences and environmental predictors was performed since SDMs must ideally restrict model calibration to accessible areas (Peterson and Soberón, 2012). Lastly, the inclusion of elevation enables model predictions not to occur in illogical areas (e.g., deep inland areas where the elevation and abiotic conditions are similar to the coastline and suitable for the species). Afterwards, the occurrence data was checked for missing (NA) values in Longitude and Latitude, and for the existence of duplicates for each species subset. Entries that matched these criteria were removed from the analysis at this point. The final dataframe was a subset for each species, and was converted into spatial polygon format for use in the SDM analysis. The post-curation number of valid entries for the analysis per species is provided in Supplementary Table 1. Lastly, a map of the global distribution of saltmarshes and mangroves ( Figure 1C of the section "Results") was also included, in order to allow the potential pinpoint of areas of pressure and potential expansion of the genus and the species with regards to the range of these two ecosystem types. This map was plotted in ArcMap v10.7.1., using two polygon layers obtained from Ocean Data Viewer 6 from (1) the Global Distribution of Saltmarshes (Mcowen et al., 2017), and (2) the World Atlas of Mangroves (Spalding et al., 2010).

Predictor Variables
For the present study, three atmospheric climate variables and one geographical variable were chosen as predictor variables. These included: (1) monthly average minimum temperature ( • C); (2) monthly average maximum temperature ( • C); (3) monthly total precipitation (mm); and finally, (4) elevation (m). All predictor variables were retrieved from the WorldClim database [version 2.1 7 (Fick and Hijmans, 2017)]. The choice of the three atmospheric climate variables was primarily based on the availability of said predictors on the WorldClim database for future periods (i.e., 2021-2040, 2041-2060, 2061-2080, and 2081-2100) as well as for historical climate data for use as a baseline (i.e., 1971-2000). WorldClim 2.1 provides downscaled data for these three variables for several Global Climate Models (GCM) of the CMIP6 (Tebaldi et al., 2021), and for the four high-priority scenarios which cover the range of possible pathways depending on socioeconomic choices. Specifically, these scenarios include: (1) SSP1-2.6 -which approximately corresponds to the previous scenario generation representative concentration pathway (RCP) 2.6, and assumes a "2 • C scenario of the sustainability" SSP1 socioeconomic family; (2) SSP2-4.5 -approximately corresponding to the RCP-4.5, and referring to the "middle of the road" socioeconomic family SSP2; (3) SSP3-7.0 -which is a medium-high reference scenario within the "region-rivalry" socioeconomic family; and (4) SSP5-8.5 -which refers to a "high reference scenario" in a high fossil-fuel development world throughout the 21st century, marking the upper edge of the SSP scenarios (Meinshausen et al., 2020). For the historical period, a total of 12 GeoTiff layers were downloaded per variable from the WorldClim database, each referring to the average monthly values between 1970 and 2000. From these, the mean of the 12 layers per environmental variable was calculated to obtain the three atmospheric climate historical predictor layers. Concerning data for the future periods, downscaled monthly climate data was retrieved from the available GCMs at the time of compilation, specifically from the ones containing all the variables, time periods, and pathways selected. The resulting GeoTiff predictor layers were restricted using a 50 km inland buffer applied to a shapefile of the world coastline, downloaded from Natural Earth Data, 8 with the intent of preventing SDMpredictions away from the coastal areas. This procedure was performed using ArcMap v10.7.1. Finally, the environmental predictor layers were stacked according to their period and SSP scenario to be used in the following SDM analysis. All variable layers collected from WorldClim and the shapefile used for geographical restriction were chosen at a spatial resolution of 10 min (i.e., ∼340 km 2 ). Layer processing was performed in R studio software v1.3.959 (R Core Team, 2020), and the script is supplied in Supplementary Material (see "SDM_predictors_Rscript").

Species Distribution Modeling
Species distribution analysis was performed using the "sdm" package (Naimi and Araújo, 2016) using the same procedures at both genus (Spartina spp.) and species level, with each species analyzed independently. First, an "sdmData" class object was created with each species, the historical climate environmental layer stack as explanatory variables, and two randomly obtained samples of the curated occurrence data for train and test data (75%). Given that the species occurrence dataframe included only presence data, an argument for background data of ten thousand (10,000) points per species using the method "gRandom" was employed, with removal of matching points, to generate pseudoabsence data (Barbet-Massin et al., 2012). This was performed since SDMs must account for absence data either with the incorporation of true data or pseudo-absence data (Peterson and Soberón, 2012); the sample size for randomly generated pseudo-absences was chosen based on Phillips and Dudík (2008). Specifically, the "gRandom" method randomly assigns absence points across the available geographical space, while the removal of matching points ensures that no pseudo-absence is placed over a known presence point. This "sdmData" object was then used for running the SDM analysis, using an ensemble of SDMs comprising generalized linear models (GLM; McCullagh and Nelder, 1989); random forest (RF; Breiman, 2001); and boosted regression tree (BRT; De'ath, 2007). In each case, four replicates with five permutations each were employed, using a train data of 70% and a test data of 30% for testing for model fit. For each model, the receiver operating characteristic curve (ROC) curves were plotted to evaluate model performance (see Supplementary Material "ROC.zip"). The mean performance (per species) using an independent test dataset is presented in Table 1. Also, for each model iteration in each analysis, the area under the curve (AUC), true skill statistic (TSS), and normalized mutual information (NMI) evaluation metrics were calculated, and each species' predictor response curves plotted (see Supplementary Material "Model_Eval.zip, " "RespCurves.zip"). From these ensembles, the historical map of habitat suitability was plotted for each species and the genus. By applying a threshold (i.e., the mean model TSS criteria of model evaluation, "max (se + sp), " the respective maps of binomial probability of occurrence (0 or 1) were obtained. By using the "predict" function (Naimi and Araújo, 2016), the fitted

Post-analysis
For each species and each SSP scenario, a set of post-analysis was performed on the resulting SDM outputs. Specifically, the change in habitat suitability (i.e., the difference between the predicted habitat suitability of a future period and the historical time) and the colonization/extinction potential (i.e., the same, but using the binomial presence/absence probability maps) were plotted. This provided an overview of the world areas in which conditions would improve or decay for each species and SSP scenario until the end of the century, while also following species potential colonization and extinction through all time points. From these maps, both the change in habitat suitability per SSP until 2100, and the change in presence and colonization/extinction potential for each SSP and future times were quantified. First, the habitat suitability net percentage change until 2100 was obtained by calculating the percentage of positive and negative cells from the habitat suitability change maps for each SSP between 2100 and historical times, for each species. These positive and negative percentages were then subtracted to calculate the net percentage change and, hence, the gain and loss of potential habitat suitability until the end of the century for each species. Second, to quantify both changes in the presence status of the genus/species, and changes in the potential colonization/extinction potential, cells of each raster map of probability of occurrence. Specifically, the frequency of cells with values above, equal, and below zero were counted in each map, and afterwards converted into percentages regarding the total number of raster cells in each map (see Supplementary Material archive "cell_count.zip"). The quantification of the potential colonization/extinction status allows the identification of areas where a certain species is projected to spread to or disappear from, as well as the total change in occurrence status globally, by calculating the net occurrence probability. These analyses were performed using the re-projection of the rasters into the Lambert Azimuthal equalarea projection, which preserves areas across the map's extent, and allows for raster cells to hold equal weight, something that the standard Mercator projection used in the rest of the study does not allow. Lastly, species overlap maps were plotted to identify areas of potential species co-occurrence and, as such, conflict and invasion/displacement potential. Specifically, for each time period and each SSP scenario, the presence/absence maps for each species were summed to obtain their respective worldwide species co-occurrence gradients. The potential change in species overlap was also plotted by subtracting each SSP scenario's time period to the historical climate species richness map. The R scripts for the post analysis are present in Supplementary Material (see "Post_analysis_Rscript").

Historical Habitat Suitability, Distribution, and Richness
The output map of habitat suitability, the geographic distribution of compiled datapoints and randomly generated pseudoabsences, and the plotted worldwide distribution of saltmarshes and mangroves are present in Figures 1A-C. Higher suitability values were observed in temperate and sub-temperate areas, partially overlapping areas of known salt marsh distribution ( Figure 1C) such as the North-eastern and South-eastern coastlines of the American continent (e.g., from the coast of Florida to Nova Scotia, and in the region north of the Amazon river's basin and the mouth of Rio de la Plata); the English Channel and the North Sea coastlines of France and the Netherlands in Europe; South-eastern Australia coastlines in the Bass Strait and the Tasman Sea, and the Northern coastlines of New Zealand; and finally in the coastline of China and Japan in the areas of the Yellow Sea, the Sea of Japan, and South China Sea. Lower values were observed to higher latitudes typically associated with polar, subpolar, and higher latitude temperate areas, as well as in the Red Sea and Gulf of Aden, and the coastlines of Peru and the Gulf of California, and other areas typically associated to mangrove habitats.
Concerning potential species co-occurrence (Figure 2), most of the areas mentioned as having higher suitability of habitat for the entire genus correspond to areas where three or more species are predicted to co-occur historically . Occurrence "hotspots" (i.e., five species, in yellow) are predicted to potentially occur along the Eastern coastline of the United States; as well as in Europe, specifically in the Gulf of Biscay and the western coastline of the Iberian Peninsula; scattered areas around the Mediterranean sea; the mouth of Rio de la Plata and the Chilean coast in South America; Northern New Zealand and the Chinese coastline in the northern region of the Yellow sea. Also, four species can potentially co-occur in several areas worldwide, including the Gulf of Mexico; the Argentinian coastline South of Rio de la Plata; the Iberian Peninsula and the English and French coastlines of the Celtic Sea; the Indian coastline of the Bay of Bengal and the Gulf of Thailand and South China Sea; the Yellow Sea and the North-eastern Japanese coastline; and also in the coastlines of Papua New Guinea and South-eastern Australia. Most areas predicted to possess four or more species overlap the known distribution of salt marshes, although the model did predict the presence of high historical species richness in areas of current mangrove distribution, between approximately 12 • N and 20 • S such as the coast between Venezuela and French Guiana in South America; the Western and Eastern African coastlines, and Indonesia, Papua New Guinea, and Malaysia.

Future Change in Habitat Suitability
The predicted potential change in habitat suitability until the year 2100 for the genus Spartina and the five species of interest is presented in Figure 3. Overall, there is a positive net change in suitability at the genus level, as well as for S. anglica and S. maritima; while this net change is negative for the remaining three species. Specifically, an increasing trend between shared socioeconomic pathway scenarios is observed for Spartina spp., which exhibits values ranging from the lowest net gain of 16.7% in SSP1-2.6, to increasingly higher values of 17.1% in SSP2-4.5, 23.7% in SSP3-7.0, and 24.9% in SSP5-8.5. Considering the five species per SSP scenario, for SSP1-2.6 the clear beneficiary is S. anglica, with an increase in suitability of just over 20%, followed by S. maritima (7.13%). On the negative end of the spectrum, S. densiflora loses the most net suitability (−14.7%), followed by S. alterniflora (−10.4%) and S. patens (just over 2%). In SSP-2-4.5, there seems to be an increase in positive net changes for all five species. Specifically, S. anglica increases up to approximately 21.8% in net suitability gain until 2100, and S. maritima registers a small increase to 4.3%. For the species on the negative end of the spectrum, their net values increase relative to SSP1-2.6, to −14.4% for S. densiflora, −8.3% for S. alterniflora, and finally −1.3% for S. patens. Under the SSP3-7.0 pathway scenario, the situation changes on several fronts. There is a decrease in the net suitability change for both species on the positive spectrum (to approximately 18.3% and under 3% for S. anglica and S. maritime, respectively), although both species continue to have a net gain in habitat suitability until the year 2100 in this scenario. Also, this scenario presumes a sharp decrease in suitability for the already negatively impacted S. patens (−9.0%) and again for S. alterniflora (−10.6%), while S. densiflora sees its best value at approximately −4%. Finally, in SSP5-8.5, S. anglica registers a net gain of approximately 21.5% in suitability, while S. maritima exhibits a value of 3.2%. On the negative end, S. densiflora exhibits a value of -5.3%, while S. alterniflora and S. patens show their lowest values comparatively to the other pathways (of −12.6 and −11.4%, respectively, in both cases). The Supplementary Material provides an animated version of habitat suitability over time for each species (see "HabitatGIFs.zip").

Future Change in Occurrence Probability
Concerning predicted occurrence distribution, determined by the percentage of positive raster cells (i.e., signifying presence) in the binomial presence and absence maps, the results follow a similar pattern as before. Specifically, most species exhibit increasing values over time that are higher with increasing severity of the SSP scenario considered. For the whole genus (Figure 4A), the increase in proportion is higher for SSP3-7.0 and SSP5-8.5, with values starting at approximately 13% for all pathways in 2021-2040 and reaching approximately 17 and 18% for SSP3-7.0 and SSP5-8.5, respectively, by the end of the century. A similar increase is observed for S. maritima and S. patens (Figures 4E,F), with both species maintaining stable values until the end of the century under SSP1-2.6, and consequently increasing with SSP scenario (e.g., an increase of approximately 3.9 and 6.7% for S. maritima and S. patens, respectively, under SSP3-7.0, and of approximately 6.3 and 8.1%, respectively, under SSP5-8.5). In S. anglica ( Figure 4B) the predicted occurrence distribution exhibits the largest increase by far with SSP, with increases of approximately 0.6 and 3.4% for SSP1-2.6 and SSP2-4.5, and of 7.6 and 10.7% for FIGURE 1 | (A) Habitat suitability map for the historical climate conditions for the genus Spartina; and (B) Distribution of presence points of Spartina spp. in geographical space (blue), plotted from the data compiled from GBIF and the surveyed literature, and the randomly generated background pseudo-absences for the whole-genus analysis (white). (C) Worldwide distribution of saltmarshes (COR) and mangroves (COR), plotted in ArcMap v10.7.1., using two polygon layers obtained from Ocean Data Viewer (https://data.unep-wcmc.org) from (1) the Global Distribution of Saltmarshes (Mcowen et al., 2017), and (2) the World Atlas of Mangroves (Spalding et al., 2010). Potential historical (1971Potential historical ( -2000 worldwide richness distribution of Spartina spp. species (i.e., S. alterniflora, S. anglica, S. densiflora, S. maritima, and S. patens) obtained from summing each species' SDM-plotted historical presence and absence predictions.

FIGURE 2 |
FIGURE 3 | Net change in habitat suitability (%) until 2100, calculated from the difference between each species' percentage gain and percentage loss in suitability between the 2081 and 2100 and the historical time periods for each CMIP6 Shared Socioeconomic pathway (i.e., . Presented values refer to each Spartina spp. species (i.e., S. alterniflora, S. anglica, S. densiflora, S. maritima, and S. patens) and for the whole genus.
Frontiers in Marine Science | www.frontiersin.org  (Figures 4C,D). Specifically, the results show predicted increases of up to 1-2% over time for either pathway in both species. Refer to the Supplementary Material for an animated version of the occurrence plots for each species (see "PresenceGIFS.zip").

Future Change in Colonization/Extinction Potential
Globally, increasing values of both colonization and extinction potential are observed over time at both genus and species level (Figure 5). Although all extinction potential values also increase with pathway severity, they are lower than those referring to potential colonization. Indeed, there is a clear tendency for increase over time and SSP for all species, with S. anglica and S. patens exhibiting the greatest increase in the proportion of raster cells for potential colonization, compared to potential extinction (Figures 5B,F, respectively). Although to a lesser extent, S. maritima, S. densiflora, and S. alterniflora also exhibit a similar pattern (Figures 5C-E), although for S. alterniflora and S. densiflora this increase is particularly small. For all species and the whole genus, the net value between the colonization and extinction potential is positive (Figure 6). In SSP1-2.6, the net value is extremely low for all species considered. Specifically, the lowest values occur for S. maritima and S. alterniflora (very close to 0% net value), with the remaining species registering values of approximately 1%, with S. anglica with the highest value. For SSP2-4.5, all species exhibit higher values comparatively, with S. anglica exhibiting an increase in approximately 3.4%, and S. densiflora now at the lower end with the lowest value. For SSP3-7.0 and SSP5-8.8, there is a further increase in all species' net change values, except for S. alterniflora. In SSP3-7.0, the value almost doubles for all species compared to SSP2-4.5, except for S. alterniflora which decreases slightly. S. anglica exhibits approximately 7.6% of net gain in colonization potential, while S. patens reaches 6.7% with S. maritima not far-below, at 4%. Finally, in SSP5-8.5 all species, except for S. alterniflora, exhibit their highest values, with S. anglica reaching almost 11% net gain, and S. patens and S. maritima nearly 8.2 and 6.3%, respectively. Despite having its highest value in this pathway as well, S. densiflora falls short of a 2% net gain. Finally, S. alterniflora exhibits its lowest value, of under 0.05% net colonization, suggesting high extinction rates for this scenario.

Future Change to Species Co-occurrence
Finally, the outputs of the species overlap analysis are provided in Figure 7. Overall, the maps show increasing levels of potential species overlap in certain regions, particularly in the higher latitudes of the Northern Hemisphere, which increase with the FIGURE 6 | Net change in occurrence increment probability (%) over time (between 2100 time and the historical period) for each Spartina spp. species (i.e., S. alterniflora, S. anglica, S. densiflora, S. maritima, and S. patens) for each CMIP6 Shared Socioeconomic pathway (i.e., SSP1-2.6, 2-4.5, 3-7.0, and 5-8.5). The proportion of colonization and extinction for each species was quantified from the change in proportion of positive (gain), negative (loss), and null (no change) raster map cells between time periods. The percentage of colonization and extinction were then subtracted to obtain the net percentage change in predicted presences. SSP in question (i.e., from SSP1-2.6 to SSP-5-8.5). Specifically, in SSP1-2.6 there appears to occur medium to long term decrease in overlap (in 2021-2040 and 2041-2060, Figure 7A) in areas such as the Eastern coast of the United States, all the way to the Atlantic coast of South America, together with losses in the Mediterranean Sea and the Iberian Peninsula and South East Asia and Oceania. In this scenario, increased species overlap is observed at a lesser, more local extent, particularly in the northern Canadian shorelines, the Baltic Sea and the North Sea in Scotland's shores, and finally in the northern Yellow Sea.
Concerning the "middle of the road" SSP2-4.5 pathway (Figure 7B), decreases are observed for smaller geographical extents, primarily in the Mediterranean Sea and southern Iberian Peninsula; while overlap increases are observed particularly in the north-eastern region of the Yellow Sea, the Baltic Sea, and the north/north-eastern region of the North American shorelines. The same pattern is present, although with progressively enhanced intensity in SSP3-7.0 and SSP5-8.5 (Figures 7C,D), with enhanced reductions of overlap in the areas of the Mediterranean and southern Iberian Peninsula shorelines and southwestern Australia coastlines, as well as in areas the areas of the eastern North American coastline and north of Rio de la Plata in South America (SSP5-8.5, Figure 7D). In these mediumhigh and high reference pathway scenarios, the increase in species overlap is also enhanced and more widespread in higher latitude areas particularly in the Northern Hemisphere (e.g., in the areas of Northern Canada and Alaska, the Baltic and North seas, the Russian coastline of the White Sea, the Yellow and Eastern China seas), but also to some extent in the Southern Hemisphere (e.g., southern Australia and New Zealand).

DISCUSSION
The present study presents a global analysis on the potential response of five species of the Spartina Genus in response to climate change until the end of the century. For the majority of the species, differential responses were observed regarding the pathway scenario in question. Despite that, however, there was a clear global trend of increasing species' occurrence until the year 2100. Also, this pattern increased with increasing severity of future scenarios suggesting general range expansion with climate change. S. anglica showed the most response of all other species, and the predicted increasing habitat suitability for all scenarios explains the increase in distribution for this species worldwide, with a very high net colonization ratio compared to the other species. Habitat suitability appears to increase with latitude over time, potentially leading to the expansion of the species' known invasive distribution in the Northern Sea and the Baltic, as well as in the Chinese eastern coastline from the Yellow Sea, where it is a known invader , potentially to the Gulf of Thailand. The present study suggests that S. anglica can potentially remain as a successful invader indifferent to the future pathway scenario of climate change, although the more severe scenarios will likely favor its expansion.
Despite being one of the more prominent cordgrass invaders (Bortolus et al., 2019), S. alterniflora showed very low levels of expansion, while exhibiting some of the worst net values for habitat suitability of all species. Notwithstanding, it followed a similar pattern in northward distribution shift as other species and remained potentially present in some of its known invasive areas. This suggests that this particular species' invasive pressure could remain, although these results seem to indicate potential vulnerability to future conditions. S. densiflora, already invasive in the Iberian Peninsula (Bortolus, 2006), exhibits the same poleward movement worldwide, particularly in Europe, suggesting that the areas of the Baltic and North seas might witness future invasion pressure from this species. However, the relatively small sample size of datapoints available for this species in areas such as the Chinese coastline could induce subsampling bias to this specific area, potentially masking effects. Finally, S. patens exhibits the second highest net gain in occurrence distribution, with increasing pathway severity, despite showing large relative losses of habitat suitability. Already an invader in Europe (Baumel et al., 2016;Sánchez et al., 2019), known also by the name of Spartina versicolor Fabre, the present study suggests that this species could not only spread its invasive distribution in the European continent (northward as the other species), it also has the potential to expand from its native areas into South America, and eventually to Asia if taken there by human action.
From an environmental management perspective, the present study identifies the SSP1-2.6 as the most beneficial. Despite changes in habitat suitability being different for all species, the net results from colonization and extinction potential worldwide present a scenario where all species maintain their current status, with relatively small increments. Also, regarding species overlap change until the century's end, despite localized decrease and gain in certain regions, this is the scenario with the lowest change, which indirectly suggests less potential for change in the current co-occurrence situations hence, less chance for newer regions of invasion-related conflicts. Exceptions are, however, the areas of the Baltic Sea and the Yellow Sea, which are presently invaded (by S. anglica and S. alterniflora, respectively), suggesting that these areas could become suitable for more species. On the other end of the spectrum, the scenario representing high fossil fueldependent development (SSP5-8.5), the potential for invasion by Spartina species increases significantly worldwide. Except for S. densiflora and S. alterniflora, this scenario leads to greater net increases in all species' potential occurrence, despite net decrease in habitat suitability in some cases (e.g., S. patens). Also, most species' distribution shows the potential to increase, particularly toward higher latitudes in the Northern Hemisphere. Accordingly, this scenario indicates greater increases in species richness northward (i.e., above 50 • N) in all northern hemisphere coastlines, while the areas of the Mediterranean and Southern Europe experience a sharp drop in richness. The Iberian Peninsula is known to harbor native S. maritima and invasive S. alterniflora, S. patens, and S. densiflora (SanLeón et al., 1999;Castillo et al., 2017;Gallego-Tévar et al., 2019), with known occurrences of S. anglica as well. In this particular scenario, S. maritima exhibits a poleward migration of its distribution, while S. alterniflora, S. densiflora, and S. patens see their potential distribution increase in this area, although with minimum overlay. These results highlight this region as highly susceptible to substituting its native cordgrass by those species that are already invading, in those scenarios where climate change mitigation measures are less efficient. Not surprisingly, this is already occurring in other areas, such as the Adriatic Sea (Wong et al., 2018). Also, in this scenario of high fossil fuel development, the region of the Yellow sea, already home to significant invasion by the smooth cordgrass (Zheng et al., 2018;Mao et al., 2019), appears particularly susceptible to invasion, since this habitat seems suitable for all the evaluated species. A peculiar result of the present study is the predicted differential range expansion/loss between colder-and warmer adapted species. Specifically, the cold-adapted species, S. anglica, is predicted to undergo an overall range expansion when compared to the rest of the species, although all species appear subjected to the same poleward pressure. This could be due to the difference in availability of longer coastlines with novel, more suitable habitats, Northward from this species' range, when compared to the rest of the analyzed species. An interesting contrast arises between the decrease in habitat suitability in species like S. patens, S. alterniflora, and S. densiflora; and the increase in probability of occurrence for the same species, although in relatively small percentages, which signify an increase in range expansion. This could be due to the fact that despite the overall suitability of habitats decreases over time, there is an increase in the latitudinal distribution of suitable shores (which is visible throrgh the expansion of the species northward), so species continue to exhibit increases in predicted occurrence (although quite limited), despite a net loss in habitat suitability.
However, some caution is advised in the interpretation of these SDMs, and these results must not be taken as absolute. First, the species presence records used were mostly obtained through GBIF, which being an open-source database is liable to sampling bias which can induce autocorrelation of the model data. Indeed, the GBIF database did not feature any points for S. alterniflora in the coastline of China, and despite the inclusion of data from articles retrieved from the available literature, nevertheless the chance for subsampling of important habitats is a potential risk. Second, since SDMs assume total occupation of climatically suitable areas (Araújo and Pearson, 2005), it is likely that there is some level of over-prediction affecting all species, culminating in predictions for larger areas. Also, the sample size of specie's occurrence points used can affect the accuracy of the predictive models (Stockwell and Peterson, 2002). In this sense, although all models included at least 300 species occurrence points -which for example is greater than the minimum number required to maximize accuracy in GLM approaches (Pearce and Ferrier, 2000) -there was a significant difference in the data size for different species (e.g., over 3000 for S. anglica while other species had less than 1000), which could lead to slightly different prediction accuracies. Another issue refers to the random assignment of pseudo-absence data across geographical space, which together with potential sampling bias (i.e., the absence of occurrence points where a species exists in reality, but has not been sampled or included in the data sources) could potentially lead to misrepresentations in the projection of geographical occupancy. Lastly, other nonclimatic variables would likely benefit the model predictions since they can help increase ecological fidelity (Austin and Van Niel, 2011). In the specific case of the Spartina species, the present study does not include potential species interaction through hybridization, although species were modeled separately to obtain a clear response to the atmospheric variables and their variation over time. Other variables salinity, immersion time, and soil composition could also be useful in this regard. Indeed, the heterogeneity of the world's coasts and the differential suitability of different types of coastline for Spartina settlement (e.g., sandy beaches, rocky shorelines, and tropical mangroves) should also be addressed, since the fraction of suitable coast significantly changes with latitude (Luijendijk et al., 2018). Indeed, the models in this study predicted and projected the presence of species of the genus Spartina in areas well known to harbor mangrove habitats, as well as some areas with relatively low saltmarsh presence. The intrusion of cordgrasses into mangrove areas is not unheard of Gao et al. (2019), and if conditions are suitable, cordgrasses have the potential to colonize mudflat areas and establish the basis for new saltmarsh habitats. The present study also highlights the potential for intrusion or colonization in these new areas, potentially establishing new competitive pressures with mangrove habitats. As such, future studies should build on the present work by testing the inclusion of these and other predictors, with the potential consequence of increasing computing requirements. Despite these limitations, however, these models provide a general indication of future trends and relative habitat suitability shifts which can be useful for environmental managers (Townhill et al., 2021).
Native saltmarsh vegetation is highly valued due to the role of its constituting species as foundation species, with substantial influence on the productivity and health of biological communities and food webs. The action of aggressive invaders when deliberately or accidentally introduced to areas outside their original range of distributions has raised significant concern from ecologists and environmental managers worldwide. The present study suggests that Spartina species can potentially benefit from climate change severity, with poleward expansions in the Northern Hemisphere being predicted for most species until 2100. Increased conflict and invasion potential could occur in Northern Europe and East Asian shorelines, in areas that already suffer from increased invasive pressure. Overall, and despite their limitations, SDMs can help establish general trends in climate change ecology and inform policymakers and environmental agents to ensure the correct management of these habitats and, ultimately, ecosystems.

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.

AUTHOR CONTRIBUTIONS
FOB, CPS, JRP, and RR performed the statistical analysis. BD, IC, and EM-N compiled and provided the occurrence data. All authors participated in the writing and revision process of the manuscript.

FUNDING
The present study was funded by Fundação para a Ciência e a Tecnologia (FCT) for funding the research via the project grant UIDB/04292/2020. This work was also funded by MAR 2020 program via the Project RESTAURA2020 (16-01-04-FMP-0014). BD and VFF were supported by investigation contracts (CEECIND/00511/2017 and DL57/2016/CP1479/CT0024). PR-S was supported by FCT through a postdoctoral grant (SFRH/BPD/95784/2013). FOB (SFRH/BD/147294/2019) and CPS (SFRH/BD/117890/2016) are currently supported by Ph.D. grants funded by community funds from FCT and the European Social Fund (ESF), through the Human Capital Operating Programme and Regional Operation Programme (Lisboa 2020). JRP is currently supported byproject ASCEND -PTDC/BIA-BMA/28609/2017 co-funded by FCT -Fundação para a Ciência e Tecnologia, I.P, Programa Operacional Regional de Lisboa, Portugal 2020 and the European Union within the project LISBOA-01-0145-FEDER-028609. The National Research Foundation of South Africa through the support of the DSI/NRF Research Chair in Shallow Water Ecosystems is thanked for funding (UID 84375) JBA.