Integrating reproductive phenology in ecological niche models changed the predicted future ranges of a marine invader

Phenology of a wide diversity of organisms has a dependency on climate, usually with reproductive periods beginning earlier in the year and lasting longer at lower latitudes. Temperature and day length are known environmental drivers of the reproductive timing of many species. Hence, reproductive phenology is sensitive to warming and is important to be considered for reliable predictions of species distributions. This is particularly relevant for rapidly spreading non‐indigenous species (NIS). In this study, we forecast the future ranges of a NIS, the seaweed Sargassum muticum, including its reproductive phenology.


| INTRODUC TI ON
New environmental conditions resulting from climate changes are expected to result in shifts of species distributions and abundances, and this raises particular concerns for expansion of non-indigenous species (NIS) (Bates et al., 2014). Climate change may benefit NIS by increasing their abundance and/or accelerating range expansions (Hellmann, Byers, Bierwagen, & Dukes, 2008;Rahel & Olden, 2008).
Such changes can potentially impact both diversity and functioning of ecosystems (Edwards & Richardson, 2004). Among the factors shaping species ranges (i.e., abiotic, biotic, dispersal strategies and capacity to adapt to new conditions; Soberon & Peterson, 2005), temperature is believed to be one of the most important (Eggert, 2012;Hoek, 1982;Lüning, 1990). In marine macroalgae (i.e., seaweeds), temperature affects individual performance (i.e., photosynthesis, growth and reproduction) and determines tolerance or survival limits (Breeman, 1990;Eggert, 2012). Quoting Breeman (1990), northern boundaries of seaweed ranges are established "by low lethal winter temperatures, or by summer temperatures too low for growth and/or reproduction," while southern boundaries are established "by high lethal summer temperatures, or by winter temperatures too high for induction of a crucial step in the life cycle." These thermal thresholds are commonly related to the tolerance of the hardiest life history stage (e.g., microthallus development), and temperature requirements for reproduction and growth (Breeman, 1990). Seaweeds are thus particularly sensitive to the cumulative effect of changing temperatures and their distributional range is expected to be strongly impacted by climate change. Warming may favour opportunistic and tolerant seaweeds by increasing their competitive ability (Dukes, 2007). Substantial changes of seaweed species ranges, including NIS, are thus expected and predicted as a response to changing seawater temperature (Assis, Serrao, Claro, Perrin, & Pearson, 2014;Breeman, 1990;Neiva et al., 2015).
In this context, ecological niche models (ENMs, i.e., species distribution models/habitat suitability models) have been increasingly used in support of invasion biology studies (Barbosa, Schneck, & Melo, 2012;Jiménez-Valverde et al., 2011). Seaweeds are no exceptions and ENMs have been used to assess the likelihood of establishment success in a novel area, identify critical routes and arrival points, and forecast range expansions and contractions under climate change (e.g., Verbruggen et al., 2013Verbruggen et al., , 2009Marcelino & Verbruggen, 2015;Chefaoui & Varela-Álvarez, 2018). Most ENMs are based on a correlative approach between species occurrence or abundance data with environmental factors and/or spatial characteristics, yielding maps of habitat suitability or probability of presence (Chefaoui, Casado-Amezúa, & Templado, 2017;Elith & Leathwick, 2009;Guisan & Thuiller, 2005). Among the challenges faced when using ENMs for studying NIS, the integration of physiological data, notably to better integrate phenology in ENMs, is critical for getting more robust predictions (Kearney & Porter, 2009;Marcelino & Verbruggen, 2015).
As highlighted above, temperature is a major driver of distributional ranges in marine seaweeds (Breeman, 1990;Eggert, 2012).
Seaweed phenology, particularly the timing and duration of the reproductive period, is strongly influenced by temperature (Breeman, 1988;Deysher, 1984). In the case of NIS, this dependency of reproductive timing on temperature might impose an additional constraint for their persistence once a new region has been colonized, which could limit the expansion of NIS in the invaded range. Day length is an additional important factor that affects reproductive timing in seaweeds (Breeman, 1993;Cunningham, Guiry, & Breeman, 1993;Voskoboinikov, Breeman, Hoek, Makarov, & Shoshina, 1996) and notably the maturation of reproductive structures in S. muticum (Strong, 2003). For instance, based on an individual-based model of population growth of the invasive kelp Undaria pinnatifida, Murphy, Johnson, and Viard (2016) showed that the fertility of gametophytes is negatively correlated with day length, suggesting that gametogenesis is mostly occurring during the shorter days of winter. Differences in recruitment patterns between the northern and southern range colonized by this NIS were also explained by different responses of the different life cycle stages to both suboptimal temperature and irradiance conditions (Murphy, Johnson, & Viard, 2017). Life cycle and reproduction of many marine species are a cumulative and combined response to temperature, day length and irradiance.
In this context, the invasive brown seaweed Sargassum muticum, native to Asia, is particularly interesting to examine. Its presence outside its native range was first reported in the 1940s, along the Pacific Coast of North America. Another major event of introduction later occurred in the 1970s along the North Atlantic coast of Europe . Among the 346 introduced seaweeds reported so far at a worldwide scale (Thomsen, Wernberg, South, & Schiel, 2016), S. muticum can be considered one of the most successful invaders. As such, it has been the focus of numerous studies. It nowadays displays a very large distribution range, spanning from Mexico to Alaska in America, and from Morocco to Norway in Europe. This pseudo-perennial species thus expanded its range over extended latitudinal ranges, in less than 40 years in Europe. In addition, diverse consequences on native species and ecosystems have been documented (e.g., Eno, Sanderson, & Conservation, 1997;Britton-Simmons, 2004;Salvaterra, Green, Crowe, & O'Gorman, 2013;DeAmicis & Foggo, 2015). For instance, S. muticum can increase biodiversity in areas where macrophytes are not present yet.
In other places, it was shown to alter native ecosystem function and structure. Considering the diverse impacts documented so far, it is important to ascertain whether the current distribution is likely to continue to expand or might contract under future climate change.
Presently, S. muticum is probably close to its distributional equilibrium with the environment since it has likely invaded a large representation of its suitable sea surface temperature in the Northern Hemisphere; thus, reducing uncertainty in niche models (Jiménez-Valverde et al., 2011). In addition, it is noteworthy that, although the information mainly comes from local studies, the species seems to display latitudinal phenological variations and acclimation to a large range of environments (see references in Engelen et al., 2015). In this study, we use ecological niche modelling to predict the future range, and thus the invasion dynamics of S. muticum, but we do so with an innovative approach by accounting for reproductive phenology. The species can reproduce only via sexual reproduction, and it is selfcompatible hermaphroditic. Fertilization is external but the zygote remains attached to the apical reproductive structures (receptacles) for 1-3 days after fertilization while initiating embryonic development. The embryos usually sink and settle immediately below or within a few meters of the maternal individual. However, a modelling approach suggested that kilometre-scale dispersal can also occur for 10%-20% of embryos under moderate current speed (Gaylord, Reed, Raimondi, Washburn, & McLean, 2002). In addition, another way to disperse is through drifting thalli. Drifting thalli cannot reattach to solid substrata or to other algae but, when carrying fertilized eggs, they can provide recruits via sexual reproduction while floating (Fletcher & Fletcher, 1975;Norton, 1977). Thus, in theory the species could colonize areas outside its reproductive window, through sexually reproducing drifting thalli, or embryos kept in suspension by turbulence. In addition, this NIS can also be dispersed by human actions, as during colonization of Europe by propagules attached to oysters. Given its high dispersal potential and semi-perennial life history (living from few years up to perhaps decades, under very favourable conditions), some populations colonizing the edge of its distribution might not be able to reproduce and release propagules.
Thus, we hypothesize that the correlative models would overestimate the future habitat for the species at the edge of its distribution, where the species could colonize and grow but never reproduce. We use reproductive windows associated with sea surface temperature and day length along the latitudinal gradient of the distribution to produce more reliable forecasts, under the assumption of conservatism of niche and reproductive physiological constraints over time. Thus, we produce a hybrid model combining species range forecasts with abiotic constraints associated with the reproductive phenology of the species. We hypothesize this approach will reduce overprediction of the future invaded range, by identifying the habitats where the NIS will not persist for generations due to reproductive limitations. To our knowledge, this study is the first to include phenological data in ENMs to improve predictions of seaweed invasion dynamics under scenarios of climate change.

| Biological data and study area
The entire known distribution range of S. muticum in the Northern Hemisphere was considered in the models, thus including occurrence data from both the native (Asia) and introduction ranges  Engelen et al., 2015). The occurrences were reviewed to correct or exclude those with referencing errors (e.g., records on land belonging to herbaria data). From the determined geographic extent, cells ranging from 0 to 30 m depth were selected, as the maximum depth at which attached S. muticum is found along the coastline.
Delimiting a depth threshold for the study area was implemented using the 30 arc-seconds resolution bathymetric data derived from the General Bathymetric Chart of the Oceans (GEBCO) website (http://www.gebco.net/). For georeferencing of S. muticum occurrence data and environmental variables, a resolution of 5 arc minutes (~9.2 km at the equator, 2.4 km at 60°N) was applied. In total, 2,587 S. muticum occurrence records were obtained, which corresponded to 1,115 presence cells (Asia: 65; N. America: 185; Europe: 865). The database of occurrence records compiled for S. muticum is available in Supporting Information Appendix S1 and Figure 1.

| Environmental data
The initial set of environmental variables for present climate conditions included the monthly average of (a) sea surface temperature (SST maximum, mean, minimum and range), and surface air temperature (SAT max., mean, min. and range), and (b) salinity. Note that SST changes have been pointed as a good indicator of climate change and its impact on phenology in marine systems (Edwards & Richardson, 2004). Surface air temperature was selected as predictor for S. muticum since the NIS can be also present on the low-intertidal shore (see e.g., Andrew & Viejo, 1998). Salinity can be also influential as the initial life stages of S. muticum can be less tolerant to brackish water (Steen, 2004). The data were derived from the Bio-ORACLE dataset (Tyberghein et al., 2012)

| Current and future species distribution models
First, we tested if the NIS was close to its distributional equilibrium with the environment in the Northern Hemisphere. We compared the range of environmental values recorded for locations of occurrence of the species within its native range (Asia), with the range of environmental values recorded in the invaded regions (Europe and North America).
To obtain the current and future climate scenarios' predictions for the distribution of S. muticum, an ensemble approach was applied. We used pooled datasets (occurrences from the native and invaded ranges) to train the models as this practice enhances the reliability of NIS predictions (Broennimann & Guisan, 2008;Chefaoui & Varela-Álvarez, 2018;Jiménez-Valverde et al., 2011).
The "biomod2" package (Thuiller, Georges, Engler, & Breiner, 2016) was used to fit six presence-absence modelling algorithms: randomForest (RF). We used an equal number of pseudo-absences extracted at random as presence cells (Lobo, 2008). All models were implemented using the "biomod2" package in R (Thuiller et al., 2016). We ran five iterations for each modelling technique, thus, 30 models for each time and scenario were generated to produce the corresponding ensembles. The data were split into F I G U R E 1 Records of Sargassum muticum included in the current analysis (magenta points) showing its distribution in the Northern Hemisphere along its native range in Asia (right) and invaded European (middle) and North American (left) ranges. Only 10% of occurrences were outside the mean sea surface temperature threshold found in the literature for the reproductive window a calibration (70%) and a validation set (30%) in each of the five iterations performed for each model. The models' performance was assessed using the true skill statistic (TSS; Allouche, Tsoar, & Kadmon, 2006), the receiver operating characteristic (ROC) curve (AUC), as well as ROC derived sensitivity and specificity (Fielding & Bell, 1997). To evaluate the predictive models, a set of thresholds was tested, and the threshold which optimized the AUC and TSS scores was chosen (Thuiller et al., 2016).
The final prediction under current conditions represented an ensemble calculated through the average of binary predictions (committee averaging ensemble), which was previously demonstrated to have the best performance in predicting coastal species distribution (Chefaoui, Assis, Duarte, & Serrão, 2016). To produce the ensemble, only the models which obtained TSS >0.8 were used. To project the future distribution of S. muticum, ensembles were projected to the RCP 2.6 and RCP 8.5 scenarios by 2050 and 2100. To assess the uncertainty of future projections, a clamping mask was generated.
The clamping mask compared the values of the variables between the training range and the future scenarios, to identify the regions where the extrapolation of models can happen. The importance of each variable was assessed using a procedure similar to "random-Forest" (Liaw & Wiener, 2002). A relative importance value from 0 to 1 (with one being of highest importance) was obtained for each variable by the correlation between the full model and a model rearranged without the examined variable using three iterations (Thuiller et al., 2016).
Our initial intention was to incorporate tolerances on both SST and day length. Thus, we estimated the mean monthly day length at the coldest and warmest months of the reproductive window found by the five studies referred above. We used the NOAA Solar Calculator (https://www.esrl.noaa.gov/gmd/grad/solcalc/) to estimate the day length during the reproductive window at the locations of the studied sites. We found that the day length reproductive window ranged from 10 hr 48 min to 15 hr 54 min (Supporting Information Appendix S2). Then, we calculated the days of the year (ranging from 0 to 365) with light within the thresholds in the entire study area (latitude 8°-70°) using the "geosphere" package in R. At the northernmost latitude (70°), we found day length within the thresholds from day 70 to 104, and from day 243 to 277. For the southernmost limit (latitude 8°), this day length threshold can be found throughout the year. Given the large range found, we concluded the day length reproductive window does not seem to be a limitation for the species persistence, thus we did not use any filter for the predictions related to day length.
For the SST window, we retrieved the minimum and maximum values of monthly mean SST of the warmest and coldest months within the S. muticum reproductive window (Supporting Information Appendix S2). Based on the data collected, the reproductive SST window of S. muticum was found to be large, ranging from 10°C to 27°C. We applied these thresholds to reclassify the ensembles of mean SST pertaining to each RCP scenario by 2050 and 2100 using the "raster" package in R. Finally, we only used the mean SST related to the reproductive window of S. muticum to filter the predictions obtained by the correlative species distribution models. We applied the true skill statistic (TSS) related threshold to transform predictions of the correlative and hybrid models into binary maps. We estimated the change in suitable habitat through the various predictions as a percentage of cells. Seas. The total area increased by 2100 under RCP 8.5 is expected to constitute 61.75%, compared to the current distribution ( Figure 4).

| Species distribution models
Overall, the area of persistence of S. muticum between the present and 2100 (RCP 8.5) was 65.14%. The persistence across time was higher in the invaded regions than in the native area.
The uncertainty in future predictions for both scenarios according to the clamping masks is presented in Figure S3.2 of the Supporting Information Appendix S3. The only areas with uncertainty due to one variable of the future scenarios presenting values outside the range of the current conditions corresponded to the eastern part of the Canadian Archipelago, located north of Hudson Bay (for both RCP scenarios), as well as the areas along the Barents and Kara Seas (for the RCP 8.5 scenario). However, these regions showing uncertainty do not coincide with the predicted areas with high probability of occurrence found by our models.

| Hybrid model
The records of the NIS match closely the mean SST upper limit corresponding to the suitable reproductive period found in the review of the literature (27°C). For the mean SST lower limit, 10% of all F I G U R E 2 Kernel density plots relating the conditions found in occurrence cells of Sargassum muticum within its native (Asia) and invaded ranges (North America and Europe) with the study area (the Northern Hemisphere). The ranges of sea surface temperature (SST) were quite similar in the native and invasive occurrences. However, S. muticum showed wider ranges of salinity in the invaded regions Study area North America Europe Asia the records were located below the threshold of 10°C: 7% of the records in Asia, 8% in Europe, and 20% in N. America (see Figure 1 and Table S3.1 in Supporting Information Appendix S3). We over-   Portugal and in Morocco (where S. muticum has been only very recently reported; Engelen et al., 2015) are also projected to retreat entirely. In addition, substantial risk of spread was also projected along the US East coast from New England to Newfoundland and along the Southern coast of Hudson Bay, a wide region where the species has not been introduced up to now.
The northward shift projected for S. muticum is in accordance with projections made for kelp forests (Assis, Araújo, & Serrão, 2018) and other fucoid algae: Fucus serratus, F. vesiculosus and Ascophyllum nodosum in the North Atlantic Jueterbock et al., 2013). These studies projected that, by 2100, these temperate species will shift northwards to Arctic shores with particularly suitable habitats being found in Spitsbergen, Greenland and Canada, while currently suitable habitats below 45°N by 2200 will become unsuitable. The northward shift in S. muticum distribution demonstrates that this invasive seaweed species may thrive well in the future.
Expansion into newly suitable habitats however requires the capacity to disperse there . Effective colonization of the North Atlantic and subarctic coastal habitats may be facilitated by the advantageous life history traits (e.g., reproduction by selfing, pseudoperennial structures) and dispersal capabilities (e.g., by drifting thalli) of S. muticum , although this may be strongly dependent on oceanographic currents (Buonomo et al., 2017). In addition, as shown by its introduction in both NE Pacific and NE Atlantic, human-mediated transport of S. muticum could easily facilitate its arrival and thus establishment in these newly suitable habitats.
Although difficult to ascertain with high certainty, it is very likely that S. muticum was introduced in the past because of aquaculture trade, in particular oyster trade . Regulations are now improved regarding aquaculture transfers but the risk of accidental introductions is not null. The new areas projected here as suitable habitats could thus become colonized quite quickly through either natural or human-mediated dispersal. In its native range, S. muticum can be found in kelp or seagrass beds together with other Sargassum species and on hard substrata . In addition to its ability to invade a great variety of unoccupied substrata, S. muticum can also invade native macroalgae assemblages and seagrass meadows if small hard substrates as pebbles, rock or oyster are available to settle on (see Engelen et al., 2015). Moreover, S. muticum is pseudoperennial, with a permanent holdfast from which new axes grow annually (Arenas & Fernández, 1998;Deysher, 1984 Rising seawater temperatures are also changing indigenous seaweed species ranges (e.g., Müller, Laepple, Bartsch, & Wiencke, 2009;Assis et al., 2014), and during such range shifts they may be more vulnerable to invasions particularly, if acting in combination with changes in other environmental parameters (e.g., Brodie et al., 2014). This may lead to the redistribution of species and differentiated populations in subarctic territories with potential evolutionary consequences (Neiva et al., 2018), as well as generate cascading and irreversible impacts on local biodiversity and ecosystems.
The reproduction-based filter had a major effect on the predicted expansion of the leading edge range, even though the number of studies that we could use were limited (5 out of 36). More data from different regions would have certainly been valuable to ascertain the variance of the reproductive window over a larger range of environments. We however do not believe that this would dramatically change our results, as the reproductive window documented so far by these studies was very large. Sargassum muticum is more uncommon in its native than in the invaded ranges, thus the kernel density plots show a higher number of occurrences at cooler temperatures in the invaded ranges ( Figure 2). However, since the lowest SST reached by the species is similar between the native and invaded ranges (Table   S3.1 in Supporting Information Appendix S3), it is likely that the species was already pre-adapted to the thermal ranges of the non-native conditions. The question of shifts in phenology and in physiology requirements following introductions is also to be further addressed (Guisan, Petitpierre, Broennimann, Daehler, & Kueffer, 2014). With the data obtained so far, which mostly rely on local studies, it is uncertain to which extent physiological requirements are stable over various regions in similar environments. In another invasive seaweed, the Pacific kelp Undaria pinnatifida, Murphy et al., (2016), Murphy et al., (2017) showed that using physiological data (e.g., growth or survival as a function of irradiance or SST) obtained in the native range perfectly predicted data observed in the field in the introduction range.
Whether such conservatism also holds in S. muticum requires further studies. Similarly, ENMs rely on the assumption of niche conservatism.
Benefitting from the introduction of S. muticum in different regions, comparative analyses of the ecosystems colonized by this seaweed would be particularly interesting to carry out to better examine putative realized niche shifts following the introduction process, as an outcome of acclimatization and adaptive evolution processes (Guisan et al., 2014;Chefaoui & Varela-Álvarez, 2018).
As Sargassum muticum is a phenotypically plastic species, new conditions generated by changing climate (i.e., increasing temperature and CO 2 levels, decreasing pH) may favour this species as it might be able to increase its competitive ability (Dukes, 2007). This may affect the entry pathways of these species, as well as its colonization, establishment, and future spread (Capdevila-Arguelles & Zilletti, 2008). In addition, intraspecific differences between invasive lineages of seaweeds in environmental requirements and niche dynamics have rarely been investigated (but see Chefaoui & Varela-Álvarez, 2018). Deeper understanding of the genetic and epigenetic factors that contribute to the invasive success of S. muticum may further improve predictions of its future spread and identify potential impacts on native ecosystems.
Up to now, projections were only made to forecast changes of S. muticum in two dimensions (i.e., horizontal/latitudinal distribution and abundance). Future research would benefit from a modelling approach that would incorporate vertical depth/intertidal height shifts, as occurring in deeper habitats or lower on the shore is a documented response of seaweed to seawater temperature changes (e.g., Pearson, Lago-Leston, & Mota, 2009); also documented for demersal fish (Dulvy et al., 2008). Range shifts and expansion in subtidal habitats could similarly occur in subtidal seaweeds. Because the upper and lower limits of seaweed distribution may be controlled by many distinct factors (e.g., Hurd, Harrison, Bischof, & Lobban, 2014), variations in other parameters as a result of climate change may result in further distributional shifts (Harley, 2011). The first consequence of long-term increases in temperature and intertidal thermal and desiccation stress (Harley, 2003), may be that the upper limits of intertidal seaweeds will shift downwards (Harley & Paine, 2009), as shown for other fucoids (Pearson et al., 2009). This may however not be the case for seaweeds that resist thermal stress by the protective effects of desiccation, for which the upper fast-drying habitat might be the most favourable (Mota, Engelen, Serrao, & Pearson, 2015).
In conclusion, this study emphasizes the importance of consid- As this species is highly tolerant and opportunistic, it may potentially displace native species populations in the regions where its distribution will be expanding. In addition, as it has been demonstrated to be a very good disperser, it may colonize these areas before range shifting natives do. In contrast, some habitats currently occupied by S. muticum (e.g., its Asian native range) will no longer be suitable and may be occupied by other seaweed species that may benefit from newly available habitats. This may result in further impacts and cascading effects on local ecosystems .

ACK N OWLED G EM ENTS
We thank Eric Treml and anonymous referees for their thoughtful comments. This study was made possible by the Erasmus Mundus