Previous studies have amplified the cox3 sequences for S. horneri populations ranging from Hokkaido, Japan (43°12′N) to Zhanjiang, China (20°53′N), including 54 sites in Japan (Uwai et al., 2009; Watanabe et al., 2019), 14 sites in Korea (Byeon et al., 2019), and 8 sites in China (Hu et al., 2011). In total, the mtDNA cox3 of 996 individuals yielded 62 haplotypes, which were downloaded from NCBI¹ (Supplementary Table S1). These sequences were aligned and trimmed using the MEGA v6.06 tool (Tamura et al., 2013). To evaluate the relationships among haplotypes, a median-joining network was built using NETWORK v4.5.1 software (Bandelt et al., 1999). Then we used the jModeltest (Posada and Crandall, 1998) to determine the most appropriate model for sequence evolution; according to the Bayesian Information Criterion (BIC), this was revealed to be Hasegawa-Kishino-Yano with a gamma distribution (HKY + G). We conducted a bootstrapped (1000 replicates) maximum-likelihood (ML) analysis using the PhyML v3.0 program (Guindon et al., 2010) and used the neighbor-joining (NJ) method in MEGA v6.06 (Tamura et al., 2013) to construct the phylogenetic trees. For these, the congeneric species Sargassum muticum (GenBank accession no. KJ938301.1), S. natans (no. KY084907.1), and S. fluitans (no. KY084909.1) were chosen to serve as the out-groups.

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⁵ genealogies and a burin-in period of 1 × 10⁶ 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).

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)², the AlgaeBase³, 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 geo-referenced 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.

Environment variables were downloaded as raster layers from Bio-ORACLE⁴ 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₂ concentration of 421 ppm by 2100, whereas the latter is a high emission “business-as-usual” scenario, with CO₂ reaching a concentration of 936 ppm by 2100 (Moss et al., 2010). For present-day time, the greenhouse gas CO₂ 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.

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 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.

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 moderate-to-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 China (=haplotype 3, in Hu et al., 2011) and South-Pacific-Japan (=haplotype D, in Watanabe et al., 2019; Figure 1B, Supplementary Table S1). A haplotype in the center of a star-like topology is often regarded as being an ancient haplotype, and H29 (AB430550) in lineage III fit this description and predominated among the samples from South-Pacific-Japan (=haplotype A, in Watanabe et al., 2019), Korea (=haplotype 1, in Byeon et al., 2019), China and Sea of Japan (=haplotype 7, in Hu et al., 2011; Figure 1B). We found the H50 (AB430580) in lineage II among samples from North-Pacific-Japan (=haplotype 39, in Hu et al., 2011), South-Pacific-Japan (=haplotype 2, in Byeon et al., 2019), and Korea (=haplotype O, in Watanabe et al., 2019; Figure 1B, Supplementary Table S1).

The divergence time between lineage I and II was c. 0.175 Mya (95% HPD 0.067–0.301 Mya) (Supplementary Figure S5), corresponding to the Saale glacial period (0.130–0.191 Mya). The estimated split time between lineage I and III was c. 0.205 Mya (95% HPD 0.035–0.414 Mya), and the divergence time between lineage II and III was c. 0.266 Mya (95% HPD 0.094–0.432 Mya) (Supplementary 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.

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).

The final predictive MAXENT model had high performance (AUC = 0.97; TSS = 0.6; Sensitivity = 0.62; Specificity = 0.98). The difference between training and testing AUC (AUC Diff.) was 0.0128 (Supplementary Table S4), without overfitting risk.

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 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).

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 (Hu et al., 2015), 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, 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-to-moderate 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.

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). Non-etheless, by 2100, the forecasted suitable areas of S. horneri differed significantly between the two carbon emission scenarios. Hence, plans for mitigating CO₂ 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).

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 warm-temperate 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).