Climate‐induced range shifts shaped the present and threaten the future genetic variability of a marine brown alga in the Northwest Pacific

Abstract Glaciation‐induced environmental changes during the last glacial maximum (LGM) have strongly influenced species' distributions and genetic diversity patterns in the northern high latitudes. However, these effects have seldom been assessed on sessile species in the Northwest Pacific. Herein, we chose the brown alga Sargassum thunbergii to test this hypothesis, by comparing present population genetic variability with inferred geographical range shifts from the LGM to the present, estimated with species distribution modelling (SDM). Projections for contrasting scenarios of future climate change were also developed to anticipate genetic diversity losses at regional scales. Results showed that S. thunbergii harbours strikingly rich genetic diversity and multiple divergent lineages in the centre‐northern range of its distribution, in contrast with a poorer genetically distinct lineage in the southern range. SDM hindcasted refugial persistence in the southern range during the LGM as well as post‐LGM expansion of 18 degrees of latitude northward. Approximate Bayesian computation (ABC) analysis further suggested that the multiple divergent lineages in the centre‐northern range limit stem from post‐LGM colonization from the southern survived lineage. This suggests divergence due to demographic bottlenecks during range expansion and massive genetic diversity loss during post‐LGM contraction in the south. The projected future range of S. thunbergii highlights the threat to unique gene pools that might be lost under global changes.


| INTRODUC TI ON
In coastal marine environments, climate oscillations during the last glacial maximum (LGM, c. 21 thousand years ago [ka]) drove coastline configurations and hence the disjunction of marine glacial refugia (Assis et al., 2017;Hu et al., 2011;Maggs et al., 2008). These refugia are isolated areas uniquely buffered from the drastic fluctuations of glacial climates and expanded to encompass not only past climate shifts, but also contemporary global warming and fluctuating sea levels (Stewart et al., 2010). The unique physical characteristics (e.g. climate dynamics and topography) surrounding climatic refugia enabled 'relict' populations to persist and retain important eco-physiological, phenological and adaptive variability (Morelli et al., 2016). However, the locations of such refugia are insufficiently understood for the northwestern Pacific coasts, where potential ancient populations persisting throughout cold and warm extreme environments may retain higher and unique genetic variability . Importantly, populations along species' range edges can have unequal adaptive capabilities to environmental shifts (King et al., 2019). Thus, quantifying evolutionary footprints of past climate-driven range shifts is key to predict the consequences of future climate scenarios, as threats to rich or endemic gene pools (Assis et al., 2017;. Coastal brown algae are ideal models to test the causal relationships between genetic diversity patterns and climate-driven range shifts . Numerous studies revealed variable roles of climate-driven range shifts, geographical isolation and genetic bottlenecks in driving phylogeographical structure and lineage differentiation in the Northeast Atlantic Neiva et al., 2015;Provan & Maggs, 2012) and Southern Hemisphere (Fraser et al., 2009). The Northwest Pacific, despite being a marine biodiversity hotspot (Briggs & Bowen, 2012) with rich seaweed diversity and endemism (Kerswell, 2006), has not yet been sufficiently studied to answer whether and how species shifted geographical ranges in response to past climate changes. Also, little information resides on the forecasted demographical consequences under future warming scenarios. It is noteworthy that the Northwest Pacific had been heavily influenced by the LGM, profoundly changing sea-level patterns and the configuration of continental shelves and margins (Wang, 1999), as well as reshaping historical demography and population connectivity (Hu et al., 2011(Hu et al., , 2015. However, no range shift estimates were yet produced, linking present-day species' genetic variation with spatial distributions across times. Sargassum thunbergii (Mertens ex Roth) Kuntze is a coldtemperate brown alga endemic to the Northwest Pacific. Historically, this species was common along the entire Pacific coast of China (Tseng, 1983), but southern populations (e.g. latitudes <21°N) have diminished over recent decades (Liu, 2017). The species can drift on the continental shelf west of the Kuroshio Current upon thallus detachment from holdfasts and potentially perform long-distance dispersal (Qi et al., 2017). Recent phylogeographical studies revealed the role of oceanic currents in driving population genetic connectivity of S. thunbergii  and identified a north to south genetic break along the coast of China, hypothetically shaped by sea-level changes during the late Quaternary . To date, no study has addressed the influence of past climate changes on the distribution of genetic pools of S. thunbergii nor estimated potential genetic losses under contrasting scenarios of future climate change.
The present study addresses two alternative process-based evolutionary hypotheses for S. thunbergii, based on both the previously reported patterns of genetic structure, and also on the sharp geographical isolation of the Northwest Pacific marginal seas during the LGM (Wang, 1999): (i) S. thunbergii persisted in demographically isolated northern and southern refugia in the Northwest Pacific during Quaternary ice ages; (ii) S. thunbergii survived in southern glacial refugia and expanded northwards following the glaciers retractions, causing homogeneous low diversity regions across the wide northern range and leaving behind most endemic diversity at rear edges. To address these hypotheses, microsatellite markers were used to detect fine-scale phylogeographical structure and latitudinal diversity patterns along the Northwest Pacific. Additionally, independent species distribution modelling (SDM) (Elith & Leathwick, 2009) and approximate Bayesian computation (ABC) were used to infer potential highlatitude range expansions and low latitude contractions, from the LGM to the present, including the warmer mid-Holocene (MH, c. 6 ka) period (Mairesse et al., 2013). Together, this information allowed testing whether past climate-driven edge-of-range shifts have influenced present genetic diversity hotspots and if these key areas are at risk of disappearing under projected global warming (Morelli et al., 2016).

| Study system
Sargassum thunbergii is found in the Northwest Pacific, as far north as Hokkaido, Japan and south to Hong Kong and Guangdong, China (Tseng, 1983), geographically spanning approximately 22° in latitude (44°-23°N; Figure 1). This species displays marked shore zonation at the intertidal to low supratidal level. It has a perennial holdfast and seasonal variations in growth of stipe, fronds, branches and vesicles, related to temperature shifts in the Northwest Pacific (Koh et al., 1993). Some Sargassaceae in China mature at 19-21°C, reaching a peak in abundance in mid-May (23.5-25°C), and perishing in late June (27.5-30°C) (Zou et al., 2006). Sexual reproduction is predominant in the population dynamics of S. thunbergii and usually needs increased resource allocation, which conversely leads to decreased allocation to vegetative regeneration from holdfasts (Chu et al., 2011).

| Sampling, genetic diversity and phylogeographical clustering
To characterize current phylogeographical diversity and genetic structuring of S. thunbergii, populations were sampled throughout its entire distribution in 35 locations between 2013 and 2016 ( Figure 1). These locations can be geographically divided into Japan-Pacific coast, the Sea of Japan, the Yellow-Bohai Sea and East China Sea, according to different water masses (Table 1). Coastal distances between sampling localities ranged from c. 2 to 2779 km. Sampling, specimen preservation and genomic DNA isolation followed Li, Hu, Gao, et al. (2017).
Eleven species-specific polymorphic microsatellite loci were used to genotype all 35 populations (717 individuals) as detailed in a previous study . Among all genotyped localities, 22 populations from China were retrieved from , and the remaining 13 populations from Korea and Japan were screened in this study. Allelic scores were checked manually for quality and consistency in STR AND (Toonen & Hughes, 2001) using the 350 ROX™ size standard (Applied Biosystems).
The number of genetic clusters was inferred using STRUCTURE 2.3.4 (Pritchard et al., 2000) without a priori population assignment and allowing admixture. For each possible number of clusters (K = 1-10), the analysis ran with 5 independent simulations using a full-length of 8×10 6 iterations and a burn-in of 10 6 iterations, which were determined to be sufficient to obtain consistent results. The estimation of K used the log probability of the data Pr(x/K) and the ΔK criterion (Evanno et al., 2005).
To compare genetic diversity and differentiation at the clustering level, ̂ , ̂ , Jost's D and F ST were also calculated in each genetic group identified by STRUCTURE analyses. F I G U R E 1 (a) Map of the study area including the latitudinal range limits of Sargassum thunbergii. (a-c) Colours depict the genetic subdivision based on structure for the two-uppermost levels of subdivisions over space (k = 2 and k = 7). Population numbers (1-35) on the map and plot structure are the same as in Table 1 1 1 1 1 1 1 K=2 K=7 K=2 These were nutrients (as nitrate and phosphate), salinity, sea ice TA B L E 1 Sites, coordinates and population sizes of Sthunbergii thunbergii analysed in this study, and standardized allelic richness (̂ ), standardized private diversity (̂ ) and gene diversity (H e ) inferred from microsatellites thickness and sea surface temperature (maximum and minimum).
Georeferenced occurrence data for the whole distribution of S.
thunbergii were compiled from the Global Biodiversity Information Facility (GBIF), the Ocean Biogeographic Information System (OBIS) and scientific journal articles (Table S1). Environmental data for the surface of the ocean were downloaded from Bio-Oracle (Assis, Tyberghein, et al., 2018). To reduce the negative effect of spatial autocorrelation (see Segurado et al., 2006), the correlation of environmental variables within the range of occurrence records was determined as a function of geographic distance with Mantel tests under 1 × 10 4 permutations (e.g. Boavida et al., 2016).
The optimal parameters of BRT models (learning rate, number of trees and tree complexity; De'ath, 2007) were tuned by implementing a cross-validation framework that partitioned occurrence records (both presences and pseudo-absences) into 10 independent latitudinal bands (Assis et al., 2017). This approach aimed to reduce the effect of overfitting and increase the potential for temporal and spatial transferability. Performance was evaluated with the area under the curve (AUC) and sensitivity (true positive rate) (Araújo, & New, 2007;Assis et al., 2017).
The potential distributional range of S. thunbergii was hindcasted (LGM, MH) and forecasted (two contrasting emission scenarios RCP26 and RCP85) with information about phylogeographical structure (i.e. within-taxon niche structure; Pearman et al., 2010). In this process, models were developed for each geographical region identified as genetically differentiated by dividing occurrence records into distinct datasets based on the identified genetic groups (e.g. Assis et al., 2016). Final predictions were reclassified to binomial responses (presence-absence maps) with a threshold maximizing specificity (true negative rate, based on pseudo-absences) and sensitivity (e.g. Assis et al., 2016Assis et al., , 2017. To assess potential niche differences between the two genetic groups, we used a probabilistic method using a Bayesian framework to determine pairwise niche overlap (see Swanson et al., 2015 for details). An additional niche similarity method using 10 4 randomizations tested whether the niche regions of the two genetic groups are more similar to each other than expected by chance (i.e. niche similarity, Broennimann et al., 2012).

| ABC inference of demographic history
The obtained phylogeographical structure and SDM results further allowed us to test a few alternative hypotheses (see Introduction) concerning the evolutionary history of S. thunbergii in the Northwest Pacific. These analyses used demographical modelling and approximate Bayesian computation (ABC) as implemented in DIYABC 2.1.0 (Cornuet et al., 2014). The first hierarchical level of structural analysis revealed a south to north genetic break (see Results); we thus conducted a relevant ABC analysis to test how these two groups established the current distribution pattern through glacial survival and postglacial colonization. We tested three broad-scale vicariance hypotheses: (i) S. thunbergii survived in a single southern refugium and expanded northward postglacially; (ii) it survived in a single northern refugium during glaciation and expanded southward postglacially; (iii) the northern and southern groups survived as separated glacial relics, respectively, from an extinct common ancestor. The divergence between the two groups could have occurred before or after the LGM, thus generating six competitive scenarios ( Figure S1). For each scenario, the current populations were set to time t 0 , and the divergence from the most recent common ancestor (MRCA) was set to t 1 (post-LGM) or t 2 (pre-LGM). The mutation rate of microsatellites (10 −4 -10 −3 ) was inferred from the phylogenetically close species Fucus (Pereyra et al., 2009). The prior values of effective population size for the southern and northern groups are listed in Table S2.
The ABC analysis generated 10 6 simulations for each scenario.
The relative posterior probability of scenarios was estimated using a weighted polytomous logistic regression method on the 1% of the simulated data sets closest to the observed data sets (Cornuet et al., 2014). Taking into account the potential lack of confidence in ABC model selection (Robert et al., 2011), we calculated the posterior predictive error over 1000 data sets to empirically evaluate the power of the model to discriminate among scenarios. In addition, we chose the mean number of alleles, mean gene diversity and mean allele size variance across loci for single sample statistics as well as the F ST values, mean index of classification and (δµ) 2 distance for two sample statistics. in the southern limits from the core-northern limit cluster (Figure 1a).

| Genetic diversity patterns and lineage structure
The second level (the optimal value of ΔK, Figure S2) divided the core-northern cluster into 6 sub-clusters of which 3 were endemic to the Yellow-Bohai Sea, and the remaining occurred along the Korean and Japanese coastlines (Figure 1b  indicated that the 7 genetic clusters were significantly differentiated from each other, particularly G2, which included the highly diverged Tateyama and Chita (Figure 2d). The curves derived from the linear regression of N a , A e , H o and H e over latitudes consistently indicated a slow but distinct decline from the northern limit (>39°N) to the centre (39°N < 28°N), and to the southern limit (<28°N) of the species range ( Figure S4). No isolation-by-distance pattern was found in the entire distribution range of S. thunbergii ( Figure S5), because P-value was statistically nonsignificant (data not shown).
Differentiation among populations within genetic clusters was greater in the range centre, followed by the northern and southern range limits, respectively (Tables S4, S5).

| Range shifts of Sargassum thunbergii from past to future
From 231 initial records of occurrence compiled from the literature and databases, a final set of 81 records resulted from considering the minimal noncorrelated distance of 15 km ( Figure S6). The BRT algorithm shows that the distribution of the species in the northern region was mostly explained by sea ice thickness (72.53%) and maximum ocean temperatures (MaxOT, 26.17%), and in the southern region by minimum ocean temperature (32.30%), followed by MaxOT (31.88%), phosphates (17.14%) and nitrates (14.48%) ( Figures S7, S8, Table S5).
The northern and southern groups in S. thunbergii have different inferred tolerance limits to salinity and nitrates ( Table 2) The posterior distribution of the probabilistic niche overlap metric between the northern and southern genetic groups is 37.3%, while the overlap between the southern and northern groups is 38.1% ( Figure S10). Niche similarity test retrieved a similar overlap of 39%, with a p-value of 0.01; rejecting the null hypothesis of climatic niches being more similar to each other than expected by chance.

| Demographic history in a model-based framework
Model-based inference at the first hierarchical level pointed to the southern group having persisted during the LGM, from which the northern group was derived through stepping stone colonization (Scenario 3; Figure S1). The southern source scenario during the LGM was supported with the maximum posterior probability (0.9696, 95% CI: 0.9642-0.9749) relative to almost zero support for the other five scenarios, including colonization from the northern source and for each deriving from an extinct common ancestor separately ( Figure S1). The estimated divergence between the southern and northern groups occurred c. 4790 generations ago (95% CI: 1400-11200), corresponding to a timeframe around 5-10 ka (95%

| Latitude-wide genetic diversity patterns in coastal marine floras
The theory predicts that leading and rear edges of a species distribution are expected to have low genetic diversity, which at the rear edge might be low within populations but with high genetic differentiation among populations due to high genetic drift, geographical isolation and selection (Eckert et al., 2008). However, this pattern was not verified in S. thunbergii. Rather, cryptic refugia, longdistance dispersal and postglacial expansion seem to have created a complex pattern of distribution range shifts and genetic diversity.
Such processes are common in seaweeds (Assis et al., 2017;Hu et al., 2011;Li, Hu, Gao, et al., 2017;Neiva et al., 2015), but were not yet understood for the complex shorelines of the Northwest Pacific.
In the Northeast Atlantic, the leading and rear edge populations of coastal floras along the same coastline exhibit contrasting genetic diversity patterns that reflect past history rather than present habitat gradients. For example, the seagrass Zostera marina is characterized by lower genetic diversity at the northern and southern range limits, and richer diversity hotspots in the centre of its distribution (Diekmann & Serrão, 2012), yet, the underlying in several other fucoid and kelp species, for which both genetic diversity within populations and differentiation among populations decrease northward from a southern range hotspot that still retains maximal diversity despite recent local bottlenecks, lower population density and higher habitat fragmentation Neiva et al., 2015). These studies stress that genetic variation patterns of coastal floras at the southern range limit can be structured by both contemporary (e.g. range dynamics and dispersal) and historical events (e.g. population bottlenecks, extinctions, colonizations and survival during ice ages) (Hampe & Petit, 2005).
Our results showed that S. thunbergii populations at the southern range limit retain less diversity ( Figure S4) and, despite being genetically isolated from more northern populations, the differentiation levels among populations within this region was the lowest (Tables S4, S5). Diversity was higher in the northern edge limit and central region, as expected considering that they comprise multiple genetic groups. Such pattern of genetic structure differs considerably from recent observations in the Northeast Atlantic Diekmann & Serrão, 2012;Neiva et al., 2015).
The unique genetic variation in S. thunbergii likely stemmed from glacial survival in the southern range, with subsequent postglacial colonization northward. This key event of the species diversification was inferred to occur in the postglacial period of 5-10 ka (ABC analysis). While significant global sea-level rise had an earlier start at ~14.5 ka (Clark et al., 2009), only at ~10 ka did the coastal morphology and the patterns of sea surface temperatures resemble those of present-day conditions (Shintani et al., 2008;Zheng et al., 2016).
In fact, the region comprising the southern range S. thunbergii had a slower warming rate during the last deglaciation, when compared to other regions of the world (Shintani et al., 2008). At the same time, selective adaptation to local environments (e.g. the northwardflowing warm Kuroshio branch current and the southward flowing cold Oyashio current along the Japan-Pacific coast) could have potentially contributed to mutation accumulation and hence fine-scale diverged lineages of S. thunbergii in the northern range.

| Glacial persistence and post-LGM range shifts shaped genetic uniqueness
The low sea levels during the LGM devastated most near-shore habitats in the Northwest Pacific (e.g. the Yellow-Bohai Sea) and imposed drastic range shifts in the extant coastal species along their geographical ranges (Wang, 1999). Generally, populations persisted locally in glacial refugia or were recently colonized from periglacial areas (Provan & Bennett, 2008). In this study, two hierarchical analyses consistently divided S. thunbergii populations south of 31°N into a unique genetic cluster (Figure 1). This cluster, along with its relatively low levels of genetic diversity at both the population and clustering levels ( Figure 2, Table 1), suggests the effect of a postglacial expansion from a presently extinct southern glacial refugium (e.g. the Hainan Island), or bottlenecks and drift owing to population reductions exposed to extreme conditions (Hu et al., 2018;Provan & Bennett, 2008).
Surprisingly, the Yellow-Bohai Sea harboured three genetic clusters (G4-G6) in the second hierarchical analysis (Figure 2c,d). The intertidal to lower supratidal distribution of S. thunbergii implies that it might not have persisted the harsher environments (e.g. ice scouring and freezing) of the Bohai Sea during the LGM, but suitable habitats might still have been available in the Yellow Sea and the west Pacific (e.g. the Ryukyu archipelago) (Figure 3a). The clusters G4-G6 likely arose from post-LGM expansions from glacial and/or periglacial ancestral relics in some areas, such as Jeju Island (Lee et al., 2012) and/or the Okinawa Trough (Hu et al., 2011(Hu et al., , 2015, driven by the Kuroshio current system . and Fucus serratus, Maggs et al., 2008). But this presumption was rejected by ABC analysis (Figure S1), which instead showed that the northern group, including G1-G3, derived from the post-LGM northward colonization of the southern group. Although the lumping of all of the northern populations together could be affecting the inference, the population size in the northern group is significantly larger than the southern (Table S1) Sulawesi, Indonesia, including the Philippine Archipelago (Figure 3a).
Such a model hindcast might have not overestimated such distribution during the LGM considering that (i) the sea surface temperatures in the South China Sea were 8-13°C during the LGM (Zhou et al., 2007) and did not exceed the current growth and survival temperature limits of S. thunbergii (4-29°C in Maizuru Bay, Japan (Umezaki, 1974); 4.2-25.8°C at Padori, Korea (Koh et al., 1993)); (ii) previous studies showed that this cold-temperate species can occur along the entire coast of China, including the southernmost Guangxi and Hainan provinces (Titlyanova et al., 2014;Tseng, 1983). During the MH, the distributional range of S. thunbergii was predicted to shift northward into the Sea of Japan and the Pacific coast of northern Japan where rising sea levels brought in the warm Kuroshio Current (Park & Park, 2015). The predicted pattern of southern decline and northern expansion of niche suitability from the LGM to the present could have also occurred in additional benthic species along the Northwest Pacific coastlines (Polovina et al., 2011).

| Conserving unique gene pools under climate warming
Genetic diversity patterns often reflect historical extinctions, dispersal and colonization events, rather than processes resulting from current ecological conditions. Inferring past distributional shifts and genetic diversity levels can provide the requisite information for contemporary conservation and management (Morelli et al., 2016), particularly for coastal floras that are vulnerable to marine climate change (Wernberg et al., 2018). In this study, the predominant role of range-wide shifts from the LGM to the present in shaping the genetic structure and diversity of S. thunbergii highlights the importance of niche conservation under climate change. This is particularly important for species like S. thunbergii that survived past climate extremes in habitats that remained suitable in the long-term, but are currently in risk of disappearing, and with it so will those unique genetic pools, as species shift distributions poleward (Lavergne et al., 2013).
The IPCC reported that the global temperature has warmed 0.85°C over the past 130 years and will be at least 1.5°C higher than in the preindustrial era by the end of 21st Century (IPCC, 2013 (Tseng, 1983). Our future predictions under warmer climate scenarios anticipate the southern range limit of S. thunbergii populations (e.g. East China Sea) even further isolated from the core-northern limit, owing to fragmented habitats and environmental pressures, particularly under the greenhouse gas emission scenario of RCP85 ( Figure 3e).
These changes caused by warming can be magnified by superimposed local stressors that might act synergistically in deteriorating the future habitat conditions for the species.
Alongside thermal conditions, the SDM models further underlined the role of nutrient and salinity regimes in explaining the southern range edges (Table S5). Continental runoffs with high nutrient concentrations from agricultural areas are already a documented concern along the coastlines of the China Sea (Breitburg et al., 2018), which can reduce habitat suitability due to local eutrophication. The additional impact of such stressors in the future distribution of S. thunbergii can be further magnified owing to the projected precipitation regimes for the region, which are suggested to increase both in frequency and intensity (Guo et al., 2018). The increase in heavy rains might further impact negatively future habitat conditions, particularly if salinity levels drop below physiological tolerance limits for any life stage, such as the most vulnerable reproductive cells.
The ongoing range shifts create multiple challenges for conservation and resource management (Lima & Wethey, 2012). The northern range shifts may first allow S. thunbergii to colonize the colder highlatitude environments to the north of Hokkaido, Japan, producing unpredictable effects on the recruitment, structure and interaction of native coastal marine species in these ecosystems (Bertness et al., 1999). Also, the unique genetic richness of S. thunbergii is under risk of erosion and loss in the northern periphery, including the Japanese interpreted the data and wrote the manuscript. All authors approved the final version of the manuscript.

Molecular data sets, genetic estimates and Species Distribution
Modelling input files are available at Dryad Digital Repository: https://doi.org/10.5061/dryad.0gb5m kkwm to be completed after manuscript is accepted for publication.