Climate-Driven Range Shifts of Brown Seaweed Sargassum horneri in the Northwest Pacific

Climate-driven shifts of coastal species’ ranges constitute a key factor shaping both the vegetation composition and biodiversity of coastal ecosystems. Although phylogeographic studies relying on genetic data have shed light on the evolutionary history of macroalgae along the Northwest Pacific (NW-Pacific) coast, their distribution dynamics are less understood and require interdisciplinary examination. Here, we used a combination of species distribution models (SDMs) and genetic data to explain the causes of population persistence and genetic differentiation and their consequences for the brown seaweed Sargassum horneri in the NW-Pacific. In this region, the phylogeographic structure of S. horneri was analyzed by screening 72 populations spawning across the entire coastal distribution and obtaining their mtDNA cox3 data. Population genetic structure and SDMs based on paleoclimatic data consistently revealed that southern coasts of the Sea of Japan, North-Pacific-Japan, and the northern part of Okinawa Trough might have served as potential refugia for S. horneri during the Last Glacial Maximum (LGM). Furthermore, we projected the distribution dynamics of S. horneri under future climate scenarios. The range of S. horneri was predicted to move northward, with a significant loss of suitable habitat, under the high emissions scenario (RCP 8.5). By contrast, projected range shifts were minimal under the low emissions scenario (RCP 2.6). Furthermore, North-Pacific-Japan was projected to be long-term persistence habitat for S. horneri under future climatic conditions, thus including this area in conservation planning could help mitigate for climate change implications. Our results enable a better understanding of the impacts of climate change on the spatio-temporal distribution of macroalgae and how this can inform coastal management and marine conservation planning.


INTRODUCTION
Understanding the impacts of global warming on seaweed communities is now a central challenge for marine species conservation. Due to anthropogenic emissions of greenhouse gases, the major consequence of climate change is an increase in both atmospheric and sea surface temperatures (SST) (Foo and Byrne, 2016), with SSTs expected to increase by 1.2-3.2 • C by the year 2100 (IPCC, 2013). The Northwest Pacific (hereon NW-Pacific) is among the most important and vulnerable marine regions for biodiversity, having lost most of the biomass because of global warming (Zhang et al., 2019;Sudo et al., 2020). Over the last three decades, from the Bering Strait to Southwestern Japan, the coastal SST has risen at a rate of 0.24 ± 0.11 • C per decade, while along the coasts of the Chinese marginal seas, the temperature has increased much faster (by 0.5 • C per decade) (Lima and Wethey, 2012), suggesting that global warming is far more complex than a simple monotonic increase of average temperatures. Identifying climate refugia and areas that harbor high biodiversity is essential to operationalize and enforce marine conservation under ongoing climate change (Rilov et al., 2019).
Canopy-forming seaweed has experienced noticeable regression because of climatic changes and anthropogenic activities (Casado-Amezúa et al., 2019). Canopy-forming species, such as large brown seaweeds support elevated levels of biodiversity through their complex structure (holdfast and surrounding substratum), thus the regression has negatively affected these associated marine species (Verdura et al., 2018;Sudo et al., 2020). Climate change has also altered the composition of canopy-forming seaweeds, with implications for ecosystem functions and services (Leclerc et al., 2015;Casado-Amezúa et al., 2019). For example, Pessarrodona et al. (2019) investigated how climatic-driven substitutions of two foundation species, Laminaria hyperborea (cold-affinity) and L. ochroleuca (warm-affinity), influenced community, trophic structure and functioning through processes associated with the cycling of organic matter (including biomass production, detritus flow, herbivory, and decomposition).
Distribution patterns of genetic diversity are often related to past demographic processes and range shifts (Dolby et al., 2016;Assis et al., 2018b). In Europe, the predicted south-tonorth recolonization-expansions in response to paleoclimatic oscillations were verified in several seaweeds, e.g., Fucus serratus (Hoarau et al., 2007), Palmaria palmata , and Ascophyllum nodosum (Olsen et al., 2010). By contrast, the postglacial expansion of warm-temperate seaweed in the NW-Pacific is affected by both complex coastal topography and ocean currents (Li et al., 2017b). Genetic diversification of marine species in this area is often related to the marginal sea basins , in that a changed sea level due to paleoclimatic oscillations has caused marginal seas to split apart (Wang, 1999). When the marginal seas reunited due to postglacial sealevel rise, populations that survived glaciation started to expand outward, with coastal currents (e.g., China Coastal Current and Kuroshio Current) possibly helping to accelerate gene flow between regions. So populations in this region have undergone multiple range contractions and expansions (Boo et al., 2019;Ng et al., 2019). Ancient isolation in separated glacial refugia is a key factor shaping the phylogeographic pattern of marine species in the NW-Pacific . Such ancient refugia were often rich in unique phylogeographic lineages that are potentially essential for their adaptation to environmental change (Gurgel et al., 2020). Most of these refugial populations, especially those at low latitude ranges, will likely become extinct under varying scenarios of climate change, resulting in the loss of this ancient genetic diversity (Assis et al., 2018a).
Climate refugia could provide comparatively stable environments wherein macroalgae could persist (Ban et al., 2016;Rilov et al., 2019), and these may be identified by species distribution models (SDMs) run under climate extremes, such as the Last Glacial Maximum (LGM, 20 Kya) and the warmer Mid-Holocene (MH, 6 Kya). SDMs are powerful tools for predicting the suitable distribution area of species and to understand how they may respond to climate change (see Renaud et al., 2019). Studies applying marine-based SDMs were predominantly carried out in the temperate North Atlantic (45%) followed by studies done at the global scale (11%) and those focused on temperate Australasia (10%) (reviewed in Robinson et al., 2017). In stark contrast, species' distribution dynamics and corresponding environmental drivers are far less understood for the NW-Pacific.
Sargassum horneri is a warm-temperate and habitat-forming species found in the coastal ecosystem of the NW-Pacific, where it provides substrate for epiphytic communities and habitat for fish and invertebrates (Komatsu et al., 2014). However, the natural resource of S. horneri is heavily threatened by human activity (e.g., over-harvesting) and climatic changes (e.g., ocean warming). As a prominent species forming the "golden tide" in the NW-Pacific, S. horneri can stay on the sea surface for ca. 1-5 months (Yatsuya, 2008). Ocean currents such as the Kuroshio Current could facilitate the long-distance dispersal of multiple marine species from south to north (Yamasaki et al., 2014;Li et al., 2017a). In particular, non-buoyant seaweed and/or invertebrates could hitchhike on floating rafts formed by positive buoyant seaweeds (e.g., Sargassum, Durvillaea, Macrocystis), to colonize to newly available suitable habitats (Fraser et al., 2013;Boo et al., 2014;Guillemin et al., 2014;van Hees et al., 2019). Thus the abundance and distribution of positive buoyant seaweeds can strongly influence the expansion of coastal species under climate change (Macreadie et al., 2011;van Hees et al., 2019).
Any range shifts of habitat-forming species found would have important consequences for community constitution and ecosystem functioning, with implications for marine management and conservation. In this study, we assessed how past climate changes produced significant biogeographical shifts and predicted the distribution of potential climatic refugia for S. horneri. We further estimated the threat posed on ancient refugia and rear edge populations by future climatic change under two contrasting RCPs (Representative Concentration Pathway scenarios).
Analysis of molecular variance (AMOVA) conducted in Arlequin 3.5 (Excoffier and Lischer, 2010) was used to assess the spatial partitioning of genetic variance among lineages. The software IMa2 (Hey, 2010) was used to estimate divergence time (t) of lineages. The final runs included a total of 1 × 10 5 genealogies and a burin-in period of 1 × 10 6 steps. The generation time, g (1 yr) and the mutation rate of cox3 (2.6-4.1%/Myr) (Uwai et al., 2009) were used to convert the output units into years. The departure from population demographic equilibrium was tested using Tajima's D (Tajima, 1989) and Fu's F S (Fu, 1997) in Arlequin 3.5 (Excoffier and Lischer, 2010) with 1,000 coalescent simulations. Under the assumption of neutrality, negative values represent population expansions, while positive values are considered a signature of recent bottlenecks (Excoffier and Lischer, 2010). Mismatch distributions implemented in Arlequin 3.5 were also used to test demographic history (Rogers and Harpending, 1992). Recent demographic expansion is characterized by smooth unimodal distribution (Rogers and Harpending, 1992).

Study Area and Occurrence Records
Our study area ranged from the Leizhou Peninsula in Guangzhou, China, to Sakhalin Bay in Russia (20-55 • N, 109-150 • E), thus covering the entire distribution of S. horneri in the NW-Pacific. We collected presence-only occurrence records of S. horneri from the Global Biodiversity Information Facility (GBIF) 2 , the AlgaeBase 3 , the journal articles, and regional experts (Supplementary Table S2). In this way, a total of 171 occurrence records of S. horneri in the NW-Pacific were finally compiled ( Figure 1A, Supplementary Table S2). Presence records were georeferenced onto a map using ArcGIS v10.2, at a resolution of 0.1 • ; any overlapping records within a cell were only considered once, so as to reduce spatial autocorrelation in the niche modeling.

Environmental Data Layers
Environment variables were downloaded as raster layers from Bio-ORACLE 4 at a 5-arcmin resolution. We selected 13 relevant variables that could affect the distribution of macroalgae in their natural environment (Supplementary Table S3). In particular, extreme environment conditions, such as the long-term average (LTA) of the hottest summer and the coldest winter month for sea surface temperature (SST) and the respective LTA of the maximum (summer) and minimum (winter) month for surface salinity (SSS), were also included. To avoid overfitting the models to the occurrence records, we used the R package 'MaxentVariableSelection' (Jueterbock, 2015) to excluding a set of environmental variables with a relative contribution score < 5% or a correlation of > 0.7 with other variables (Supplementary Figures S1, S2, Table S4). Collinearity of the retained environmental variables were assessed with the variance inflation factor (VIF < 10). VIF was assessed using the R package usdm (Naimi et al., 2014).
Environmental variables (temperature and salinity) for the LGM were derived from PaleoMARSPEC data layers (Sbrocco and Barber, 2013). To project these into future climate conditions, we chose two widely recognized representative concentration pathways (RCPs): RCP 2.6 and RCP 8.5, the former being the most optimistic climate change scenario with a predicted CO 2 concentration of 421 ppm by 2100, whereas the latter is a high emission "business-as-usual" scenario, with CO 2 reaching a concentration of 936 ppm by 2100 (Moss et al., 2010). For present-day time, the greenhouse gas CO 2 concentration used was 410 ppm (d.d. May 2018; www.esrl.noaa.gov) (Westmeijer et al., 2019). Apart from temperature, the monthly maximum and minimum average sea surface RCP scenarios also included simulations of salinity changes arising from alterations to the hydrological cycle. Both SST and SSS were each derived for the mid-century (2040-2050) and end-century (2090-2100) in the two RCP scenarios (from Bio-ORACLE).
All the environmental layers were gridded to a 1-arcmin resolution (approximately 0.016 • at the equator), by using bilinear interpolation. In all cases, the projected distribution area was restricted in space by a 17.5-km buffer zone from the coastline.

Species Distribution Models (SDMs)
To generate SDMs for S. horneri, we applied Maximum Entropy Modeling in MAXENT v3.3.3 (Phillips et al., 2006) and Genetic Algorithm for Rule-set Prediction (GARP) model in openModeller v1.2 (Muñoz et al., 2011). These algorithms are compatible with the low numbers of available information on the natural occurrence of our study species. We randomly selected 20% of the distribution records for training, reserving the remaining 80% for testing. MAXENT models the environmental conditions in the distributional range of S. horneri according to 10 000 background sites positioned randomly along the coastline. Then MAXENT was run 15 times, and each run was convergent within 1000 iterations. According to the output results, the variables with a relative contribution score less than 5% in the model were excluded. For the GARP model, we ran it with 10 000 replicate enforcing a low omission (E = 5%), and selected 10 best models out 100 runs. Default settings were otherwise used for all other parameters in each model. A jackknife test was used to rank the importance of each environmental variable in both models; the variable with the smallest contribution was dropped, and this process repeated until only one variable remained which we inferred was most important variable for predicting suitable habitat. The final model was run using all occurrence data and selected variables to project the suitable habitat within the study area.
The models' algorithm performance was evaluated using the predicted area under the curve (AUC), provided by the ROC curve (Receiver Operator Characteristic) (Phillips et al., 2006), ROC-derived sensitivity (presences correctly predicted) and specificity (absences correctly predicted) (Fielding and Bell, 1997), and the true skill statistic (TSS) (Allouche et al., 2006). For Frontiers in Marine Science | www.frontiersin.org AUC, models with value > 0.90 were considered "excellent" and in the range 0.7-0.9 "reasonable predictions". For TSS, models with value > 0.8 were considered "excellent" and in the range 0.4-0.8 "good" (de la Hoz et al., 2019). MAXENT showed better performance (AUC = 0.97; TSS = 0.60) in the validation than did GARP (AUC = 0.87; TSS = 0.70), so we used MAXENT to generate the distribution models using all presence points of the seaweed species. Then we projected these models onto ancient and future climatic layers. Additionally, model overfitting was calculated as the difference between training and testing AUC, and higher values (AUC Diff.) indicate over-parameterizing (Jiménez-Alfaro et al., 2018).
Distribution maps with threshold values represent habitat suitability (range: 0-1) were built up by MAXENT. To illustrate the changes to suitable habitat in response to climate changes, we compared the future suitable area and southern ranges with its present distribution (Figure 2). The probability grid data were projected in UTM (Universal Transverse Mercator Grid System) to calculate the area of suitable habitat (species' presence probability > 0.5) location arising under different periods and RCP scenarios.

Genetic Structure
The haplotype network revealed deep genetic splits among populations along the NW-Pacific (Figures 1A,B). The ML and NJ phylogeographic trees separated the 62 haplotypes into three genetic lineages (labeled here I, II, and III), with moderateto-high support for each node as indicated by their bootstrap values (Supplementary Figure S3). Lineage I dominated the populations from the North-Pacific-Japan; Lineage II dominated coasts of Pacific-Japan and Korea; Lineage III occurred in the populations from South-Pacific-Japan, China, Korea, and Sea of Japan (Supplementary Table S1, Figures 1A,B). In particular, 98.5% samples from South-Pacific-Japan grouped in lineage III, and the rest samples were found in lineage II. Notably, all the haplotypes found in China and the Sea of Japan belonged to lineage III ( Figure 1B). For AMOVA analysis, 84.10% of genetic variance occurred between lineages (F ST = 0.84170, P < 0.0001; Supplementary Table S5).
Most haplotypes were represented by individuals located in one biogeographic region, apart from the three haplotypes H5, H29, and H50 ( Figure 1B, Supplementary Table S1). Of the latter, H5 (JF461023) in lineage III was found in samples from  Figure S5). For lineage I and II, Tajima's D and Fu's Fs values were statistically non-significant, and the mismatch distribution of the two lineages showed no signs of recent demographic or range expansion (Supplementary Figure S6). The Tajima's D and Fu's Fs values estimated for lineage III were significantly negative (p < 0.05) (Supplementary Figure S6). In addition, unimodal distribution was observed in the mismatch distribution for lineages III (Supplementary Figure S6), indicating a demographic expansion scenario.

Environment Variables and Model Performance
The model with the lowest AIC was built based on the following three uncorrelated environmental variables (Supplementary Figure S1, S2, Table S4): mean primary productivity, LTA of maximum sea surface temperature (MaxSST) and LTA of minimum surface sea salinity (MinSSS) ( Table 1). For GARP, variables with an omission error above the 95% basic omission error (=6.33) could change the potential distribution area significantly ( Table 1). The ranking of relative contribution of the three variables was the same in both MAXENT and GARP models ( Table 1). Mean primary productivity was the most important variable (44.7% model contribution, Table 1) in discriminating suitable from non-suitable habitats. MaxSST had a contribution of 28.2%, and MinSSS had a contribution of 27.1% in Maxent (Table 1). Response curves showed that mean primary productivity was positively correlated with the habitat suitability of S. horneri whereas MinSSS was negatively correlated with it (Supplementary Figure S4). The MaxSST in the range of 24.5 • C-28 • C was generally predicted as optimal for the species (Supplementary Figure S4).

Ancient, Present, and Future Projections of Habitat Suitability
When projected back in time, to the conditions of the LGM (20 Kya), both models indicated the southern Sea of Japan were potential refugia for S. horneri (Figure 1C). In addition, the northern Okinawa Trough and North-Pacific-Japan were suboptimal (with probability of 0.5). During the MH (6 Kya) period, S. horneri began to expand outside refugia as sea levels rose, moving westward to the Bohai Bay and northward to the coasts of Hokkaido, Japan, the Russian coast of the Japan Sea, the Sakhalin Bay and also around the Kuril Islands ( Figure 1D). The present-day model's results were consistent with the species' current distribution (Figure 2). These projections showed that habitat suitability was high along the coast of East China Sea and Bohai Bay, the western and southern coasts of the Korean Peninsula, the western part of Kyushu Island and the coasts of Honshu Island (Figure 2). Compared with the Pacific coasts of Japan, a lower coverage of suitable area was detected along the Sea of Japan. Since we did not include seashore substrata in the seaweed's niche modeling, some of the predicted suitable habitat was inconsistent with occurrence records reported for some areas. For example, in our modeling, high habitat suitability was detected along the coasts of Jiangsu MaxSST, long-term average of maximum surface temperature; MinSSS, long-term average of minimum surface salinity.
Province of China (30 • -35 • N), but this area lacks the rocky substrata that are essential for macroalgae establishment and growth. At the present-day time, S. horneri was projected to have 16.75% suitable habitat (Figure 3). During the LGM, the suitable habitat was contracted to 15.17% of studied area (Figure 3), ranging from c. 32 • -41 • N (Figure 4). By 2050, under the stabilizing scenario RCP 2.6, S. horneri will lose 13.43% of its current suitable habitat, but by 2100, suitable area will arise to 17.14% (Figures 2, 3). Notably, the suitable habitat of S. horneri increased along the coasts of Korea under RCP 2.6 (Figure 2), whereas that sustained at its southern boundaries changed little (Figures 3, 4). Under RCP 8.5, the more severe climate change scenario, by 2050 the highest habitat suitability was found along the coasts of Central-Pacific-Japan, and the southern and eastern coasts of Korea (Figure 2). Projections also indicated that parts of the coasts of China had suitable habitat for S. horneri. By 2100, however, S. horneri populations should be able to colonize new suitable habitat available in the north yet they will lose all their suitable habitat along the Chinese and Korean coasts, with only the coasts of Central-Pacific-Japan and northern Sea of Japan having the highest habitat suitability (Figure 2). The suitable habitat of S. horneri was contracted to11.87% and 2.08% of studied area by 2050 and 2100, respectively (Figure 3). The northward shift of southern boundaries was predicted to be ca. 1581.89 km (ca. 14.26 • ) by 2100 (Figure 4).

DISCUSSION
Sargassum horneri, a warm-temperate brown algal species, forms ecosystem-structuring underwater forests. Genetic divergence was detected in S. horneri populations along the coasts of NW-Pacific, suggesting multiple ancient refugia (e.g., North-Pacific-Japan, South Sea of Japan, and northern Okinawa Trough). These refugia were supported by ancient distribution of S. horneri during the LGM. Furthermore, the loss of habitat at the southern boundaries was predicted by 2100 under high CO 2 emission scenario (RCP 8.5).

Ancient Distribution and Lineage Divergence of the Seaweed
Within a species, lineage divergence reflects the occurrence of an allopatric distribution (Ni et al., 2014), and paleoclimatic oscillations are regarded as being one of the most crucial factors shaping the phylogeograhic structure of coastal species (Zhang et al., 2019). Highly structured phylogeographic patterns were detected in most marine species that occur along the NW-Pacific coasts, such as the seaweed Saccharina japonica (Zhang et al., 2019), Chondrus ocellatus , the yellowfin goby Acanthogobius flavimanus (Hirase et al., 2020), and the barnacle Chthamalus moro (Wu et al., 2015). Sargassum horneri populations in the NW-Pacific were also characterized by strong genetic variation among regions (Hu et al., 2011). In our study, the divergence time between lineages ranged from 0.175 Mya to 0.266 Mya, much earlier than the LGM (0.018-0.020 Mya).
Populations from North-and Central-Pacific-Japan grouped in lineage I and II, leading to a deep genetic split from other populations (Figure 1). Hu et al. (2011) presumed the existence of the northern cryptic refugia and our projections for the LGM support this view, in that North-Pacific-Japan (ca. 35-40 • N) may have served as a glacial refugium for S. horneri. Similarly, an ecological niche model also predicted that another seaweed, S. japonica, survived the LGM in this same area (Zhang et al., 2019). Lineage II also occurred along the coasts of Korea. The Kuroshio Current, bifurcating at the southwestern Japan, acts as a hydrological barrier between populations in eastern Korea and Pacific coastline of Japan; thus, the occurrence of lineage II in both regions likely resulted from a much earlier mixture (e.g., pre-LGM) rather than from recent genetic exchange.
Based on our projections, the southern Sea of Japan represented a suitable distribution area during the LGM (20 Kya). In addition, S. horneri in South-Pacific-Japan and Korea was highly structured based on cox3 sequences, implying long-term persistence in this region. During the LGM, the SST was about 3 • C-5 • C lower than it currently is in Japanese waters (Oba et al., 1991). Since S. horneri can grow rapidly under low temperatures (10 • C -15 • C) (Yoshida et al., 1998(Yoshida et al., , 2004, the lower SST in the past might have only shortened its growth period or slowed its growth rate instead of causing its large-scale extinction. Although the northern Okinawa Trough was projected to be low-tomoderate suitable area for S. horneri (Figure 1C), we assumed a high probability of ancient populations once existing in this area. During the LGM, the Yellow-Bohai Sea and East China Sea were reduced to Okinawa Trough, and Seto Inland Sea was terrestrial (Wang, 1999). The private haplotypes along the Chinese coasts are strongly indicative of refugia in the Okinawa Trough, as the elapsed years since the LGM (ca. 20 Kya) are too short for the formation of such a large amount of genetic variation. Furthermore, it is evident that the populations in lineage III have experienced a population expansion (Supplementary Figure S6). In line with our interpretation, Cheang et al. (2008) inferred that the "Paleo-coast" of the East China Sea could have served as glacial refugia for a congener species Sargassum hemiphyllum, leading to the vicariance of its ancestral populations.

Range Shifts and Environmental Drivers
Small environmental changes can be magnified at the rear edge of species' range shift, which could adversely affect their distribution and diversity (Assis et al., 2016). In our study, the projected losses of suitable habitats in southern boundaries (namely South-Pacific-Japan and the East China Sea) for S. horneri were detected under RCP8.5. Similarly, southern boundaries of S. muticum were predicted to shift northward, from the coast of China and Korea to the Russian coast of the Sea of Japan, the Sakhalin Bay, and around the Kuril Islands (Chefaoui et al., 2019). Further, Westmeijer et al. (2019) estimated that some macroalgae species (Fucus serratus, Laminaria digitata, Gracilaria gracilis, Furcellaria lumbricalis, etc.) would shift 110 to 635 km from the coast of the Iberian Peninsula toward northern Brittany by 2100. The dispersal and distribution patterns of marine species are also driven by ocean currents (Li et al., 2017a). Future climate-driven oceanographic changes are likely to strengthen the Kuroshio Current, which may increase the rates of northward expansion of marine species in the NW-Pacific (Wilson et al., 2016). Nonetheless, by 2100, the forecasted suitable areas of S. horneri differed significantly between the two carbon emission scenarios. Hence, plans for mitigating CO 2 emissions could be adopted to prevent the local extirpations and range shifts of this and similar species due to ocean warming.
Our findings showed that primary productivity was the most important driver shaping the distribution of S. horneri. Primary productivity, itself closely associated with nutrient concentrations, available light, and water transparency, could reliably discriminate suitable from non-suitable habitat for S. horneri populations. In particular, the concentration of nutrients seems to have been a relevant factor driving the population dynamics of S. horneri (Yoshida et al., 2001;Qi et al., 2017). Response curves showed that mean primary productivity was positively correlated with the habitat suitability of S. horneri. This is consistent with primary productivity projected as being a major variable shaping the current distribution of seaweeds Undaria pinnatifida (Báez et al., 2010), Macrocystis and Durvillaea (Bernardes Batista et al., 2018).
Extreme temperatures are likely to impact the survival of a critical phase in the life-cycle of marine species, which is more important for determining their distribution range (Martínez et al., 2018;Figueroa et al., 2019). Summer and winter SST variables are often related to shifts of southern and northern range boundaries, respectively (Neiva et al., 2015). In our study, summer temperature (MaxSST) contributed more than did winter temperature (MinSST) when projecting the distribution pattern of S. horneri. Most temperate seaweeds at the rear edge often suffer from summer isotherms above their maximum thermal thresholds (Eggert, 2012), causing their local extirpations. Summer temperature was deemed the dominant factor affecting the projected distribution boundaries of most canopy-forming seaweed in the Northwest Atlantic (Wilson et al., 2019), northern Japan (Sudo et al., 2020), and Australia (Martínez et al., 2018). The thermal thresholds of species are related to their life stages, thus integrating phenological trends should improve our ability to predict future distributional shifts at both regional and global scales. For example, Chefaoui and Serrao (2017) demonstrated that the total range area of S. muticum was expected to increase by just 1.63% when the reproductive temperature window was explicitly considered in model projections.
Seasonal salinity variation was also detected along the coasts of China due to the changes in oceanic water budget and circulation caused by river discharge (Yellow River and Yangtze River). For S. horneri, its sporophytes recruit in winter and grow rapidly from November to March (Xu et al., 2016;Zhang et al., 2018). MinSSS was projected to be a critical driver determining the distribution range of this species, perhaps by influencing its embryos development and germlings' growth. In two earlier studies, the adults of S. muticum tolerated a wider range of salinity than did its germlings (Steen, 2004), which might also apply to S. horneri. This low salinity also affects the survival and growth of S. hemiphyllum, and it neither found in the Bohai Sea nor the coast north of the Yangtze River (Cheang et al., 2008).

Consequences of Range Shifts and Implications
Previous studies have noticed that canopy-forming marine species are suffering particularly strong declines worldwide (Assis et al., 2017;Wilson et al., 2019). According to work by Krumhansl et al. (2014), 38% of the world's kelp forests have vanished in the past five decades. Sargassum horneri is one of the most important habitat-forming species in the NW-Pacific, where it annually forms large-scale floating mats (Li et al., 2020). These floating islands of biomass accumulate at the continental shelf and the Kuroshio Front, providing both shelter and food for juveniles of many fish and invertebrates (Komatsu et al., 2014;Yamasaki et al., 2014). Based on satellite data, large patches of floating S. horneri were first found near Zhejiang Province waters, and then transported northeast (Komatsu et al., 2007;Qi et al., 2017). Due to their role in supporting biodiversity and food webs, the changes in the distribution range of S. horneri may substantially impact the species diversity of coastal and pelagic ecosystems (Yamasaki et al., 2014).
Our study combined species distribution model and genetic analyses to give insight into impacts of climate change on the spatio-temporal distribution of S. horneri in the NW-Pacific. Based on our projections, the North-Pacific-Japan region may be a suitable habitat for S. horneri under ongoing global warming, which is also supported by projections done for S. japonica (Zhang et al., 2019). The cold Oyashio current in this region may mitigate the rate of ocean warming and facilitate the growth of cold-water adapted species kelps (Saccharina, Alaria, Costaria) (Sudo et al., 2020) and warmtemperate species, such as Undaria pinnatifida (Sato et al., 2016) and Sargassum thunbergii (Li et al., 2017a). Under the RCP 8.5 scenario, our projections demonstrated a significant range contraction of S. horneri with the loss of its southern limits. Rear-edge populations may embody unique adaptive strategies, which could maintain population stability so these populations should be of special interest to conservation planners (Rilov et al., 2019).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
J-JL conceived the study and coordinated all steps in finalizing the manuscript. J-JL and S-HH carried out the data analysis and wrote the manuscript. Z-YL provided guidance on the modeling. Z-YL and Y-XB interpreted the results and contributed to the writing. All the authors contributed to both the data collection and the manuscript's revision and also contributed to the article and approved the submitted version.