Ecological niche dynamics of three invasive marine species under the conservatism and shift niche hypotheses

Marine bioinvasions are one of the main threats to biodiversity. According to assumptions based on ecological niche models, the conservatism and equilibrium of species with the environment are vital to understanding the bioinvasion process. However, these assumptions have been evaluated primarily for terrestrial species, with few examples in marine environments. We tested the niche conservatism and niche shift hypotheses in native and invaded environments and evaluated the niche dynamics and invasion stages on three invasive marine species: the algae Asparagopsis armata and Codium fragile , and the ascidian Asterocarpa humilis . We applied the identity and background similarity tests to assess the conservatism, the principal component analysis to evaluate the niche dynamics, the Gallien et al. approach to evaluate the invasion stages, and an ensemble of model to estimate potential distribution. Findings showed that the niche equivalence hypothesis was not rejected for any of the species, indicating equivalent ecological niches. Niche similarity demonstrated that niches in native and invaded ranges were not similar as expected by chance for A. humilis and C. fragile . However, for A. armata , the populations in the native and invaded areas had a very similar environmental niche. In addition, high niche stability is evident in the niche dynamics, and so as the stabilization phase of the invasion phases of the three species; thus, studying the three species supported the hypothesis of niche conservatism. These results indicated that all three species have dispersed and are in biogeographic equilibrium within their invaded regions.


Introduction
Non-indigenous species (NIS) are one of the main drivers of change in marine biodiversity and a significant biosecurity problem worldwide (Simberloff et al. 2013;Katsanevakis et al. 2014;Ojaveer et al. 2018). Consequently, studying the native and invaded geographical distributions of NIS and their invasion dynamics at different ecological and temporal Citation: Rivera R, Pinochet J, Brante A (2022) Ecological niche dynamics of three invasive marine species under the conservatism and shift niche hypotheses.
Aquatic Invasions 17 (in press) scales is crucial to generate prevention and mitigation plans and management policies for NIS (Pullin 2002;Villero et al. 2017).
Reliable biological and environmental information is essential to evaluate hypotheses and generate robust predictive models in the invasion phenomena (Ferrier et al. 2002). The current availability of global environmental databases (e.g., MARSPEC -Sbrocco and Barber 2013;Bio-ORACLE -Assis et al. 2018), geo-referenced records of species (e.g., GBIF, OBIS), and the development of the geographic information system (GIS) have allowed significant advances to study and test hypotheses related to the spatial distribution of species. Among these tools, ecological niche models (ENMs) that use spatially explicit information have been widely implemented to study biological invasions (Peterson 2003;Jiménez-Valverde et al. 2011;Jeschke and Strayer 2008).
Under the framework of biological invasions, the fundamental assumption of ENMs is that species retain their native ecological niches in the invaded area, implying their inability to adapt to new environmental conditions (niche conservatism hypothesis; Peterson et al. 1999;Peterson 2011). Thus, ENMs implicitly assume that NIS retain their niches (Peterson and Vieglais 2001;Peterson 2003Peterson , 2011Václavík and Meentemeyer 2012). However, contrary to expectations, recent evidence has suggested that some taxa do not exhibit niche conservatism but rather "niche shift" (Broennimann et al. 2007;Petitpierre et al. 2012;Pili et al. 2020). This assumption would produce a climatic mismatch between the native and the invaded areas due to a rapid adaptive response of NIS caused by selection pressures in the new environment (Broennimann et al. 2007;Gallagher et al. 2010;Oliveira et al. 2018). Likewise, the biological factors caused changes in the fundamental niche (species evolve in a new range); changes in the realized niche; disequilibrium in the native range; different biotic interactions; climate unavailability in the introduced range; or modeling issues (Zachariah Atwater and Barney 2021). In this way, short-term changes may occur in the ecological niches of NIS, which, in some cases have proven to be substantially different from their native ranges Gallagher et al. 2010;Atwater et al. 2018). These changes suggest that the new environmental conditions faced by NIS, in many cases, are not an obstacle for their establishment (Gallagher et al. 2010).
To study the niche dynamics, Petitpierre et al. (2012) developed a framework showing the three components of the niche: stability, expansion, and unfilling. Here, the stability component is the proportion of the total niche shared by the native and the invaded niches (i.e., climatic match). On the other hand, the expansion corresponds to a fraction of the invaded niche not shared with the native niche; this area corresponds to newly occupied environments in the invaded range that have not been observed in the native range . Finally, the unfilling component shows environmental conditions from the native range unoccupied by populations in the non-native range (i.e., climatic mismatch).
Another approach allows studying the phase in which a NIS is a framework developed by Gallien et al. (2012). This approach compares a global niche model considering all available occurrences (native and nonnative) with a regional niche model (i.e., occurrences in the non-native range). The occurrence can be classified into four invasion stages: stabilizing population (both models return high suitability), this means that those individuals have established or will most likely become established, sink population (both models predict low suitability), regional colonization (the global model predicts high suitability, while the regional model predicts low suitability), and regional adaptation (the regional model predicts high suitability even though the global model predicts low suitability). An essential feature of this approach is the interpretation of the proportion of stabilizing populations as an estimate for invasion success.
Contrary to terrestrial systems (e.g., Oliveira et al. 2018;Pili et al. 2020), few studies have been carried out in the marine realm to evaluate invasion processes considering the scenarios of niche conservatism and niche shift (but see Parravicini et al. 2015;Carlos-Júnior et al. 2015;Battini et al. 2019;Melo-Merino et al. 2020). Studies incorporating these scenarios could better understand the spatial patterns and ecological processes that operate in the phenomena of invasions in marine systems, which differ significantly from terrestrial ecosystems in terms of ecological components and evolutionary history. This work involved the evaluation of the invasion dynamics of three model marine species, Asparagopsis armata (Harvey, 1855), Codium fragile (Hariot, 1889) ssp. tomentosoides, and Asterocarpa humilis (Heller, 1878), and testing the niche conservatism and niche shift hypotheses. The selection of these species was due to the extensive biological and ecological information available in the literature. The species were selected using two criteria. First, the species poses a significant risk to biodiversity and, possibly, the economy; the second criterion is an extensive geographical description of the species since the models must be robust.
The red alga Asparagopsis armata is native to Australia and New Zealand (Andreakis et al. 2004). This species is considered a pest in the Mediterranean, where it has significantly impacted native communities, outcompeting key native species and becoming an economically harmful species (Kraan and Barrington 2005). Its primary dispersal mechanism is the rafting of tetrasporophyte and gametophyte fronds (Martins et al. 2019). Codium fragile is a green alga native to Southeast Asia and an invader on the coasts of the Mediterranean (Provan et al. 2005), North Pacific (Provan et al. 2008), andChile (Villaseñor-Parada et al. 2018). This alga has a series of traits that seem to favor its ability to invade new habitats, such as a high tolerance to fluctuations in temperature, salinity, light, and nutrients, as well as its ability to reproduce both sexually via gamete fusion and asexually through parthenogenesis and fragmentation (Chapman 1999;Trowbridge 1998;Mathieson et al. 2003;Harris and Jones 2005). In a recent study comparing the traits of 113 macroalgae species introduced in Europe, C. fragile was considered the riskiest species in terms of dispersal capacity, probability of establishment, and ecological impacts on the host communities (Nyberg and Wallentinus 2005). Asterocarpa humilis is an ascidian, initially identified in New Zealand and later reported in different coastal regions around the world, such as the SE Pacific (Chile), SE Atlantic (South Africa), and recently in the Northeast Atlantic (France and Great Britain; Kott 1985;Clarke and Castilla 2000;Bishop et al. 2013;Turon et al. 2016;Pinochet et al. 2017). The main dispersal mechanisms of this species are maritime transport on ship hulls and aquaculture (Clarke and Castilla 2000;Bishop et al. 2013;Pinochet et al. 2017).
We tested the hypotheses of niche conservatism and niche shift with these three model species. We also investigated niche dynamics (i.e., stability, unfilled, and expansion), we evaluated invasion stages (i.e., adaptation, stabilizing population, sink population, and colonization), geographic extension, and the main environmental drivers that explain the invaded range of these species using multiple models within an ensemble forecasting framework.

Study area
Provided that the extension of the geographic area can influence the performance of ENMs (Barve et al. 2011), the spatial scale of the study area must consider the dispersal capacity of the target species . As a proxy of each species' dispersion potential and estimation of its geographical range, we generated a buffer equal to half the mean distance between occurrences (Marshall and Strine 2019). This distance was calculated using a Behrman projected coordinate system. For Asparagopsis armata, the average distance was 13941 m, Asterocarpa humilis, 820 m, and Codium fragile 54677 m (see Supplementary material Figures S1, S2), as a proxy of their dispersion potential ). We considered the study areas and calibrated the ENMs. These areas were considered as the accessible world (M) according to Soberón and Peterson (2005). Native and invaded areas for A. armata. A. humilis and C. fragile are available in Figures S1, S2.

Occurrences
Presence records of the three model NIS were compiled from the open databases Ocean Biogeographic Information System (OBIS, http://www.iobis. org) and GBIF (Global Biodiversity Information Facility, https://www.gbif.org), accessing directly the portals cited above, as well as records available in specialized literature. Access to all the mentioned online databases occurred on December 15, 2020 (Table S1). All the macroalgae records were also cross-checked using http://algaebase.org. The records were manually reviewed in ArcGIS 10.5, excluding duplicate data, points positioned on land areas or without coordinates, and those with low coordinate accuracy (i.e., degrees only). After data filtering, 1828 records for Asparagopsis armata, 229 for Asterocarpa humilis, and 999 for Codium fragile were used for posterior analyses. In addition, to avoid spatial biases in the sampling effort, which are common when literature and databases are used (García-Roselló et al. 2015), we carried out a spatial thinning approach to eliminate records with a minimum distance of 9.2 kilometers to assure one point of occurrence per cell grid of the environmental layers. By doing so, we made sure to avoid redundancy in our analysis. Reducing the thinning distance to a value lower than 9.5 would considerably reduce the number of occurrences.
For this reason, we opted to find a balance between the number of unique, valid occurrences that were not autocorrelated. The analysis was carried out using the spThin R package (Aiello-Lammens et al. 2015). After spatial thinning, 547 records were recovered for A. armata, 92 for A. humilis, and 395 for C. fragile ( Figure 1).

Environmental database
Environmental predictors were obtained from Bio-Oracle v 2.1 (Assis et al. 2018) at five (5) arcminutes (~ 9.2 km) resolution. This database has been developed for ecological niche modeling, and they have been widely used to predict potential distributions of non-native marine species in coastal areas (Zhang et al. 2020;Derviche et al. 2021;López-Farrán et al. 2021). This dataset included 11 environmental variables (Table S2) restricted to the top surface layer and described monthly averages for the period 2000-2014. However, considering that the correlation between variables may affect the performance of ENMs (De Marco and Nóbrega 2018;Feng et al. 2019), we performed GLM with a binomial distribution to exploratory analyzes for each species, the explanatory variables were standardized (difference from the mean divided by the corresponding standard deviation). With these analyses, we evaluated the statistical significance of the 11 environmental predictors on the presence of each species. After a staggered process, only the significant variables were retained (p < 0.05) (Table S3). We then ran GLM models and got a reduced model with only significant variables. From a total of 11 environmental variables, a subset was selected, given their low value of VIF (< 3) (Zuur et al. 2010) and a correlation coefficient | < 0.7 | (see Figure S3). Likewise, following the proposal of Cobos et al. (2019), we did not choose to use heuristic methods since the performance decrease in models created with variables selected in this way (Cobos et al. 2019). Thus, the variables used for posterior analyses were the mean values of calcite, dissolved oxygen concentration, pH, primary productivity, salinity, silicate concentration, and temperature. The VIF analysis was performed using the USDM package (Naimi et al. 2014).

Niche conservatism and niche dynamics
To test the niche conservatism hypothesis and niche shift, we first calculated niche overlap, which is the area shared by the invaded and native niches. For this calculation, we estimated Schoener's D overlap index, which ranges from 0 (no overlap) to 1 (complete overlap) (Warren et al. 2008). Schoener's D calculates the suitable range for a given species based on probability distributions for inhabiting particular regions (cells), calculating niche overlap based upon species abundance in those locations (Warren et al. 2008). D indexes were evaluated under two alternative scenarios. (1) niche equivalence, that is, if the observed niche overlap is indistinguishable when the grouped occurrences of both niches are randomly reassigned, and (2) niche similarity, which estimates if the occupied niche in one geographic range (native area) is more similar to the occupied niche in the other geographic range (invaded area) than randomly generated niches. Niche equivalence was evaluated by simulating two niches based on occurrence records, randomly permuting the native and invaded niches. The Schoener's D (Schoener 1970) overlap index was estimated using 999 randomizations to generate a null distribution of randomly permuted niche overlap values (Warren et al. 2008). Niche conservatism was inferred if significant values for the niche equivalency test indicate that niches are more similar than would be expected by chance. Furthermore, significant values for the niche similarity test (in any comparison direction: native vs. non-native and non-native vs. native) indicate that niches are more similar than would be expected by chance (Warren et al. 2008). Analyses were carried out with the "ecospat" package (Di Cola et al. 2017) implemented in the R software (R Core Team 2020).
The native and invaded niches dynamics were evaluated using a PCA-env analysis . This method uses the first two axes of a PCA to define the environmental space, and the presence records are transformed into occurrence densities using a Kernel function . A total of 10,000 pseudo-absences were generated to estimate the density of available environments in the environmental space. The occurrence density and available environments densities were used to estimate an occupancy index to compare the densities of occurrence delimiting the occupied niche in the native and invaded ranges . It is important to note that the framework presented by Petitpierre et al. (2012) and the niche conservatism and shift hypotheses must be evaluated in environmental analogs . We evaluated this assumption (described extensively in the Ecological Niche Modeling section) by performing the ENMs only in geographic areas corresponding to environmental analogs (i.e., similar areas in environmental terms).

Invasion stage
To evaluate the invasion stage of each NIS, we followed the proposal of Gallien et al. (2012) that considered four categories. First, the presence of the species was predicted by the global ENMs (occurrences in the native and invaded ranges) and the regional ENMs (i.e., occurrences in the native range only). In this case, species were considered to be in equilibrium with the new environment. Second, the occurrence of the species were not predicted by the global or the regional ENMs. Under this scenario, the species occupied unsuitable habitats. Third, the occurrence of the species were predicted only by the global ENMs, but not by the regional ENMs, indicating that they were at a colonization stage. Fourth, the occurrence was predicted by the regional ENMs, but not by the global ENMs; suggesting expansion of the species into new habitats. The evaluation of the invasion stage was performed using the occurrences predicted by both global and regional models, which were correlated with the actual occurrences recorded in the invaded range (Gallien et al. 2012;Langdon et al. 2019). The analyses were carried out in the "raster" package (Hijmans 2020) with the R software (R Core Team 2020).

Ecological niche modeling
Calibration of the models To evaluate the geographic extent of Asparagopsis armata, Asterocarpa humilis, and Codium fragile, we modeled the environmental niche for the three species using nine algorithms consisting of different approaches to ecological niche modeling. We opted for ensemble modeling since this approach considers that each algorithm has a true "signal" as well as some "noise" created by error and uncertainty of the data and the structure of the model (Araújo and New 2007). Provided that the uncertainty in the distributions predicted by the different models can skew results (Araújo and New 2007), we explored a range of predictions through different ENMs to generate a consensus projection, a method known as ensemble species distribution models (ESDMs) (Araújo and New 2007). We used a weighted average as a consensus method. The utilized models were: 1) generalized additive model (GAM), 2) generalized linear model (GLM), 3) multivariate adaptive regression splines (MARS), 4) boosted regression trees (BRT), 5) classification and regression trees (CART), 6) random forest, (RF), 7) Maxent, 8) artificial neural network (ANN), and 9) support vector machine (SVM). The predictor variables used for all models were the seven (i.e., calcite, dissolved oxygen concentration, pH, primary productivity, salinity, silicate concentration, and temperature) selected to evaluate the niche conservatism hypothesis.
Additionally, we generated uncertainty maps that represent the variance among methods. Three approaches were used to generate the ENMs: 1) the models were calibrated in the native range and projected towards the invaded areas (ENM inv ), 2) the model was calibrated in invaded areas and projected into the native areas (ENM nat ); and 3) models used records of native and invaded areas together (ENM global ). All analyses were performed using the "SSDM" package (Schmitt et al. 2017). We used the area under the ROC curve (AUC) average using threefold cross-validation (two times in the training set and once the evaluation set), and repeated the process ten (10) times. For selecting the best-fitted model and for finally building the ensemble probability surface, the AUC threshold was kept at < 0.8. Scores ranged between 0.5 and 1, with those close to 1 indicating that the model has excellent performance, whereas values around 0.5 denote models that are no better than random. Following recommendations from previous studies (Barbet-Massin et al. 2012), we chose 10.000 random points restricted the background (pseudo-absences) to ecologically plausible areas. To create the final ensemble model, we calculated the average of all algorithms whose AUC scores were over 0.8. We projected the final ensemble models for each species onto a predicted distribution map. Additionally, we created uncertainty maps using the predicted probability of occurrence variance among the algorithms, where areas of agreement are characterized by the low variance between algorithms (Schmitt et al. 2017).

Performance evaluation of the models
We evaluated the performance of the ENMs using a multi-metric approach to determine the variability among estimates. Due to the criticisms received by the AUC (Fielding and Bell 1997) (see Lobo et al. 2008;Peterson et al. 2008), we additionally calculated the AUC ratio of the partial ROC index (pROC) (Peterson et al. 2008) and Boyce's index (Hirzel et al. 2006) for presence-only records. The pROC was calculated using the "ntbox" package (Osorio-Olvera et al. 2020), and the statistical significance was assessed by bootstrapping and comparing to a null hypothesis (i.e., H 0 = difference between model-predicted AUC and random AUC is ≤ 0). We used 30% of the presence records randomly chosen for resampling with 500 iterations. Significance was evaluated using the calculated AUC values and the values of pseudo-replicates following the proposal of Peterson et al. (2008). The Boyce index was calculated as the relationship between the density of occurrence records that fell within the predicted habitat thresholds and the total number of sites where the species were found and recorded (Boyce et al. 2002;Hirzel et al. 2006). Index values ranged between −1 and 1, with negative values indicating that the model predicts inadequate areas where the presence point is found, and positive values designating that the model predicts suitable areas where the presence points are located; that is, they indicate close agreement between the prediction of the model and the true or known distribution of the species (Boyce et al. 2002;Hirzel et al. 2006). The Boyce index calculations were performed using the "ecospat" package (Di Cola et al. 2017).

Assessment of model assumptions (Extrapolation risk)
Considering that the evaluation of niche conservatism/shift under the Petitpierre et al. (2012) approach must be carried out in environmental analogs ), we evaluated the risk of extrapolation. This phenomenon occurs when models generated in one geographic area and projected in another may contain new climatic values not available in the projected areas (Owens et al. 2013). There are two kinds of new conditions: 1) for an individual variable, the values may be outside the range sampled during training (i.e., strict extrapolation); and 2) portions of the environmental space can be located within the range of individual variables, but represent a new combination of predictors (Owens et al. 2013). We identified analogous environments by comparing native and invaded ranges using three approaches to assess the risk of extrapolation: a) Analysis of Multivariate Environmental Similarity Surface (MESS; Elith et al. 2010; Figure S4); this metric provides an index of environmental similarity between each pixel and the median of the most dissimilar variable (Owens et al. 2013). b) Mobility-Oriented Parity (MOP), which identifies strict extrapolation areas and calculates the environmental similarity between the calibrated and projected areas (Owens et al. 2013; Figure S4). c) Extrapolation detection (ExDet; Mesgaran et al. 2014), which identifies similar or novel environments between native and invaded areas (Mesgaran et al. 2014). Unlike the other approaches, ExDet identifies the most influential environmental variables leading to non-analogous environments between the compared areas (Mesgaran et al. 2014). NT1 ranges from infinite negative values to zero, where zero indicates no extrapolation beyond the univariate coverage of reference data ( Figure S5). The MESS and MOP analyses were performed in the "ntbox" package (Osorio-Olvera et al. 2020), and ExDet was calculated in the software ExDet v1.1 (Mesgaran et al. 2014).

Niche conservatism and niche dynamics
The first two PCA axes compare the climatic niche between native and invaded ranges, explaining the 41.85% and 18.87% of the variance in Asparagopsis armata, 43.7% and 19.85% in Asterocarpa humilis, and 50.78% and 17.46% in Codium fragile. The ecological niche overlap (Schoener's index) between the native and invaded ranges was < 0.35 for the three species (Table 1).
The niche equivalence hypothesis was not rejected (Identity test; p > 0.05) for any of the three species, indicating equivalent ecological niches between their native and invaded ranges (Table 1, Figures S6, S7, and S8). The hypothesis of niche similarity (or the background test) indicated a similarity more significant than that expected by chance for A. armata (Similarity Test; p < 0.05), indicating that the populations in the native and invaded areas have a highly similar environmental niche (Table 1). For A. humilis and C. fragile, the test indicated that niches in the native and invaded ranges were not more similar than expected by chance (Similarity Test; p > 0.05), suggesting differences in their ecological niches ( Figures S6, S7 and S8). According to the niche dynamics analysis (PCA-env approach), the three species presented high stability (A. armata 98%, A. humilis 72%, and C. fragile 50%), with similar environmental conditions occurring in both the native and invaded ranges (50%; Figure 2; Table 1). Regarding the expansion phase, C. fragile presented a 50% expansion of its niche, followed by A. humilis with 28%, and close to zero for A. armata (0.02%) ( Figure 2; Table 1). The results indicated low unfilling values for the three species (< 0.1%), revealing a low proportion of environmental space present only in the native range, which translates into few highly suitable sites in the invaded range compared to the native range ( Figure 2; Table 1).

Invasion stage
The invasion stage for Asparagopsis armata revealed that 69.6% of the sites were in the stabilization stage, 24% showed a colonization state, 4.3% represented adaptations to new environmental conditions, and 2.1% presented characteristics of sink populations (Gallien A). For Asterocarpa humilis, 60.3% of the sites were in the stabilization stage, 30.6% represented adaptations to new environmental conditions, 8.2% corresponded to sink populations, and 0.9% were in the colonization stage ( Figure 3B). Lastly, for Codium fragile, 84.6% of the sites represented stable areas, 8.4% of the populations were in the colonization stage, 6.1% corresponded to sink populations, and 0.9% were adapting to new environments ( Figure 3C).

Ecological niche modeling
The predicted geographic distributions indicated good performance mainly for the global models (i.e., calibrated with records of the native and invaded area; Table 2) concerning the reciprocal models (i.e., calibration in native area and projection in the invaded area and vice versa). For Codium fragile, high values of the Boyce index (0.96) and pROC (1.94) were observed (Table 2). Asparagopsis armata presented a better performance in the AUC (0.94) and pROC metrics (1.93; Table 3). In the case of Asterocarpa humilis, the global model presented the best performance (Boyce index = 0.994), and the reciprocal comparisons (invaded to native) presented the lowest values, with an AUC = 0.822, Boyce index = −0.860, and a non-significant pROC (Table 2). Regarding the contribution of the variables that explained Figure 2. Niche dynamics of a) Asparagopsis armata, b) Asterocarpa humilis, and c) Codium fragile. Niche comparisons were carried out between the native and non-native ranges of each species. Colored areas represent stability (blue), unfilling (green), and expansion (red). Solid and dashed lines delimit 100% and 75% of the available background environment considered in the analyses. The red arrow is the distance between centroids.   the geographical distribution of the three NIS, the results indicated that the highest contributions were provided by calcite, followed by primary productivity and salinity (Table 3), for global models as well as reciprocal comparisons ( Figures S9, S10). The geographic predictions (ensemble forecasting) indicated that for A. armata, there are areas of high suitability along the coasts of Chile, Peru, South Africa, Great Britain, the Sea of Japan, and the southeastern coasts of Australia and New Zealand ( Figure 4A). For A. humilis, the areas of highest suitability are limited to Chile, Peru, southern Australia, and New Zealand in the South Pacific, the coasts of Argentina and South Africa in the South Atlantic, Great Britain, and the China Sea ( Figure 4B). Geographic predictions for C. fragile revealed areas of high suitability in Peru, Chile, New Zealand and the Tasman Sea in the South Pacific, Argentina and Namibia in the South Atlantic, the Pacific and Atlantic coasts of the United States, northwestern Europe, the Black Sea, the Sea of Japan, and the East China Sea ( Figure 4C).
Since different algorithms can provide different predictions even when using the same dataset, the use of model ensembles allowed us to identify the uncertainty of the predictions. Indeed, the greatest uncertainty in A. armata, A. humilis, and C. fragile predictions was found in areas, far from the coast, where the observed presence of NIS has been confirmed ( Figure S11).

Climatic analogy
The MESS analysis revealed climatic analogs between native and invaded regions for the three NIS ( Figure S4). In the case of mobility-oriented parity (MOP), discrete zones of strict extrapolation were observed, that is, dissimilarity of environments between the calibration (M) and projection. For Asparagopsis armata, they appeared mainly in northern Europe, the Mediterranean, and North Africa. For Asterocarpa. humilis, few areas of strict extrapolation appeared; whereas, for Codium fragile, a more significant number of areas with environmental non-analogs were found ( Figure S4). A complementary analysis to detect non-climatic analogs (ExDet) was performed, which revealed that for the three NIS, there was no extreme extrapolation ( Figure S5). These results indicate that the niche dynamics tests (the framework proposed by Petitpierre et al. (2012)) were sufficiently robust because they met the condition of being evaluated exclusively in analogous environments.

Discussion
The literature on biological invasions indicates a controversy in that NIS predominantly present ecological niche conservatism, that is, introduced species colonize environments similar to those of their native ranges (Araújo and Pearson 2005;Wiens and Graham 2005;Wiens et al. 2010;Guisan et al. 2014) or NIS experience a niche shift in the invasion process (Broennimann et al. 2007;Petitpierre et al. 2012;Early and Sax 2014;Pili et al. 2020). Nevertheless, our work revealed that the three species, Asparagopsis armata, Asterocarpa humilis, and Codium fragile, exhibited ecological niche conservatism. Three lines of evidence support this finding. First, all three species presented a niche equivalence, the most used proxy to support ecological niche conservatism as it indicates no significant alteration of their environmental niche during invasion (see Broennimann et al. 2012). Secondly, the dynamics of the ecological niches of all three NIS revealed a high percentage of stabilization, i.e. these species presented a high shared proportion of their native and invaded niches, a reflection of niche conservatism Langdon et al. 2019). Finally, the invasion stages revealed that the stabilization phase was predominant in all three species; indicating that all three NIS are in equilibrium in their introduced ranges. The similarity test identified a significant niche similarity between the native and invaded ranges for A. armata, thus providing further evidence of niche conservatism for this species, but was not conclusive for A. humilis or C. fragile.
For Asparagopsis armata, both similarity and niche equivalence results indicated that this species occupies similar niches in its native and invaded ranges (see Table 1). It proved to be in a stage of stability in its invasion dynamics ( Figure 2A) and a stabilization phase in its invasion stage ( Figure 3A). These results could be explained by the decrease in the number of occurrences reported by GBIF and OBIS in recent years ( Figure S12), indicative of a potential decrease in the arrival of this species to new locations. Despite this, the invasion success of A. armata could be explained by its physiological advantages, such as high tolerance to low temperatures (Maggs and Stegenga 1999), which, combined with favorable physical-environmental factors, could determine the establishment of this species in diverse invaded habitats. Indeed, Martins et al. (2019) indicated that interactions with conspecific and macrophytes could determine the distribution and abundance patterns of A. armata.
Similar results were obtained for Codium fragile, showing that this species also tends to occupy similar environments in its native and invaded ranges, thus indicating a highly conserved ecological niche. Regarding its ecological niche dynamics, C. fragile has been proved to be in a stage of stability, presenting a shared proportion between its native and invaded ranges, indicating niche conservatism Langdon et al. 2019; Figure 2C). The above is reinforced by the fact that this NIS is at the stabilization phase, that is, the species previously colonized new sites and is currently in equilibrium within these new environments ( Figure 3C). Biogeographically, C. fragile has a wide worldwide distribution, being considered a highly invasive species (Roman and Darling 2007;Therriault and Herborg 2008;Rius et al. 2014). This extensive geographical distribution could be explained by its wide physiological tolerance (Madariaga et al. 2014), an important characteristic for the success of the invasion process of the littoral biota. In addition, it inhabits intertidal and subtidal environments and artificial structures (e.g., dock, aquaculture implements; Neill et al. 2006) that can favor its distribution success. Codium fragile has been recognized as the main invasive macroalgae in these systems (Castilla et al. 2005), negatively affecting the yield in macroalgae cultivation (Neill et al. 2006). However, at a local level, biotic interactions, such as herbivory and possibly competitive exclusion, could slow down the rapid invasion process, which could explain the geographic stability described in the literature (Gulliksen and Skaeveland 1973;Yamaguchi 1975;Petersen and Svane 1995), as well as that reported in this study. This stability pattern could also be explained by the decrease in the number of geo-referenced occurrences in recent years ( Figure S12), reflecting a decreasing trend in its geographic expansion. However, one attribute of this NIS is its sexual and asexual reproductive mechanisms, bestowing it with a high capacity for dispersal and invasion (Madariaga et al. 2014), indicating that this species could still colonize new environments and potentially settle in new areas.
Similar environmental niches in the native and invaded ranges of Asterocarpa humilis were also found, showing high stability in terms of its niche dynamics ( Figure 2B). A high percentage of its populations also proved to be in the stabilization stage, suggesting that it is in equilibrium in the new invaded environments. However, the ecological niche modeling with the best fit of the regional model (i.e., invaded range) with respect to the global model (i.e., native + invaded; Figure 3B) suggests that this species could present local adaptations, which could indicate it is still in late expansion. According to Pinochet et al. (2017), A. humilis is still being transported over long distances by adhering to ship hulls, allowing this species to disperse to new geographical areas where it has already been described as invasive, for example, Chile (Pinochet et al. 2017), France and Great Britain (Bishop et al. 2013). In general, both the presence of aquaculture and boats allow this species to colonize new sites even though these are not environmentally suitable (see Figure 4), representing a potential late expansion.
The ENMs revealed robust results for the three species according to the standard AUC, partial ROC, and Boyce index tests with the geographical predictions of spatial dispersion. Indeed the global models (those that use the occurrences of the native and invaded ranges) were superior in performance compared to reciprocal or regional comparisons indicating that the latter is insufficient to predict the invasion potential of a NIS (Sales et al. 2021). The superior performance of global models have previously been described in the literature, emphasizing the need to integrate native and invaded ranges occurrences to capture more significant environmental variability where NIS occurs (e.g., Escobar et al. 2018;Langdon et al. 2019;Pili et al. 2020). Therefore, we demonstrated that the use of ENMs provides a valuable tool to detect areas at risk of introduction and dispersal of NIS at a macro-spatial scale, being an efficient invasibility proxy and permitting the detection of areas that could be of interest to prevent or detect early introductions. NIS monitoring in already colonized areas is also essential to prevent future spread scenarios.
The use of multi-approach tools in this work provides an integrative view of the invasion process in three NIS of high ecological and economic importance, considering the negative impacts they cause at an ecological and economic level. It should be noted that the study of niche dynamics in marine NIS, as opposed to terrestrial species, has been poorly documented (e.g., Petitpierre et al. 2012;Goncalves et al. 2014;Oliveira et al. 2018;Herrando-Moraira et al. 2019;Hu et al. 2016), and mainly limited to the Mediterranean (Parravicini et al. 2015;Poursanidis et al. 2020). Among future research to be developed in the prediction of invasive species, the Bayesian spatial method should be used, allowing the incorporation of spatial correlation of the variables and the uncertainty of the parameters in the modeling process, resulting in a better quantification of the uncertainty and accurate predictions (Ellison 2004;Pennino et al. 2014). Bayesian spatial models may also aid data analyses with geographically varying levels of survey effort, which reduces its influence on estimates of the effects of environmental variables (Gelfand et al. 2006).
Finally, in this study we show that these three species of importance as NIS in Chile, have dispersed and are in biogeographic equilibrium within their invaded regions. Which is unlikely to have a large dispersion greater than what has already been recorded.

Supplementary material
The following supplementary material is available for this article: Table S1. Occurrence points used in ecological niche modeling. Table S2. Environmental variables used to develop ecological niche models. Table S3. Summary of the generalized linear models (GLM) for the relationship between NIS and different environmental variables. Figure S1. Study area (native) used in the ecological niche model of Asparagopsis armata, Asterocarpa humilis and Codium fragile. Figure S2. Areas invaded by Asparagopsis armata, Asterocarpa humilis and Codium fragile. Figure S3. Spearman's correlation for environmental variables used in the models of Aparagopsis armata, Asterocarpa humilis and Codium fragile. Figure S4. Results of the multivariate environmental similarity surface (MESS) and mobility-oriented parity (MOP). Figure S5. ExDet map of the extrapolated projection areas (non-analog areas). Figure S6. Statistical tests for niche comparison between native and non-native ranges for Asparagopsis armata. Figure S7. Statistical tests for niche comparison between native and non-native ranges for Asterocarpa humilis. Figure S8. Statistical tests for niche comparison between native and non-native ranges for Codium fragile. Figure S9. Native model projected onto the invaded range. Figure S10. Invaded model projected into the native range. Figure S11. Uncertainty map based on the between-algorithm variance. Figure S12. Frequency of occurrences per year (OBIS+GBIF+literature).