Agarophyton vermiculophyllum populations were collected during 2014–2018 from Hokkaido, Japan (44.05°N) to Beihai, China (21.45°N). In total, 546 specimens were sampled from 39 sites (14 from Japan, nine from Korea, and 16 from China) (Supplementary Table S1). At each site, 5–33 individuals were collected randomly with an interval distance >10 m to minimize the chance of collecting genetically identical samples. All samples were stored in silica gel. Due to the scarcity of Korean samples, molecular data of eight additional sites from Korea (62 specimens) and one site from Russia (three specimens) were retrieved from Kim et al. (2010) [marked with asterisk (^∗) in Supplementary Table S1]. In total, we analyzed 611 samples of A. vermiculophyllum from 48 populations.

Genomic DNA was extracted from 5–20-mg dried samples using the Onestep-Lysis^TM Plant DNA Isolation Kit (PosPolysaccharides & Polyphenolics-rich, Nobelab, Beijing, China) following the manufacturer’s instructions. We amplified partial mitochondrial DNA (mtDNA) cox1 region with primers 43F (5′-TCA ACA AAT CAT AAA GAT ATT GGW ACT-3′) and 1549R (5′-AGG CAT TTC TTC AAA NGT ATG ATA-3′) (Geraldino et al., 2006). PCR amplification was performed in 50 μl of reaction volume, containing 25 μl 2 × EasyTaq Master Mix (Cwbio, Beijing, China), 1 μl of each forward and reverse primer (10 μmol L^–1), 1 μl DNA template (15 ng μl^–1) and 22 μl RNase-free water. Amplification was run in a PCR Thermal Cycler (LongGene, China) with an initial denaturation at 94°C for 10 min followed by 35 cycles of denaturation at 94°C for 30 s, annealing at 50°C for 30 s and extension at 72°C for 2 min, with a final extension at 72°C for 10 min (Yang et al., 2008). PCR products were purified using a Gel Extraction Kit (Omega, Norcross, GA, United States), and bidirectional sequencing was performed on the ABI3730 sequencing platform (Applied Biosystems, Foster City, CA, United States) using another primer set 622F (5′-CCT GTN TTA GCA GGW GCT ATT ACA ATG C-3′) and 880R (5′-ACA GTA TAC ATA TGA TGN GCT CAA AC-3′) (Yang et al., 2008).

BioEdit v.7.2.5 (Hall, 1999) was used to align newly obtained sequences with previously published eight mtDNA cox1 (GenBank: EF434932–EF434339) (Yang et al., 2008). In total, 619 mtDNA cox1 sequences including eight reference sequences were imported into MEGA v.7.0.21 (Kumar et al., 2016) for a comprehensive alignment and manual editing. Four genetic diversity indices, the number of haplotype (Nh), the number of polymorphic sites (S), haplotype diversity (Hd) and nucleotide diversity (π), and neutrality tests (Tajima’s D and Fu’s Fs) were estimated in Arlequin v.3.5.1.3 (Excoffier et al., 2005). The distribution of characteristic nucleotide mutations in each cox1 haplotype was calculated in MEGA.

Bayesian inference (BI) and maximum likelihood (ML) trees were reconstructed using BEAST v.1.8.0 (Drummond et al., 2012) and PhyML v.3.1 (Guindon and Gascuel, 2003), respectively. The best-fit nucleotide substitution model of mtDNA cox1 (GTR + G, G = 0.1490) was determined based on the Akaike Information Criterion in jModelTest v.0.1.1 (Posada, 2008). For BI analysis, a Yule process tree prior with the UPGMA tree was constructed in BEAST. An uncorrected lognormal relaxed clock with default estimated rate of 1.0 was set. BEAST was run with 5 × 10⁷ Markov Chain Monte Carlo (MCMC), sampling every 1,000 steps under the GTR + G model. The MCMC output was the tree data set and checked for convergence [effective sample sizes (ESS) >200] in Tracer v.1.7.1 (Rambaut et al., 2018). To obtain the condensed tree, the MCMC output was summarized in TreeAnnotator v.1.7.5 (Drummond et al., 2012), the first 10% was discarded as “burn-in”. The ML tree was constructed based on nearest neighbor interchange (NNI). Phylogenetic support for each node was assessed by bootstrapping with heuristic analysis of 1,000 replicates. For both the BI and ML analysis, the other two members of the genus Agarophyton were chosen as outgroup (Gurgel et al., 2018), i.e., A. tenuistipitatum (C. F. Chang et B.- M. Xia) Gurgel, J. N. Norris et Fredericq (GenBank: EF434924), and A. chilense (C. J. Bird, McLachlan et E. C. Oliveira) Gurgel, J. N.Norris et Fredericq (GenBank: EF434915).

Genetic clustering analysis was inferred based on an admixture model using rmaverick v.1.0.5 (Verity, 2019) in the R package. ClustalX v.1.83 (Larkin et al., 2007), xmfaconvertor, and xmfa2struct were used to convert cox1 sequence files to the rmaverick input .str files. The values of K were set from 1 to 10, and 10 replicates were run for each K using 10⁴ burn-in iterations and 2,000 sampling iterations.

Genealogic relationships of mtDNA cox1 haplotypes were inferred using the median-joining algorithm implemented in Network v.5.0.1.1 (Bandelt et al., 1999). Intraspecific relationships of cox1 haplotypes were also assessed by minimum spanning network (MSN) in Arlequin.

Analyses of molecular variance (AMOVA) were conducted in Arlequin to detect spatial partitioning of genetic variance among sampling sites and regions. We used four criteria for population grouping: (i) Country, populations were divided into three groups: Japan, Korea, and China, and the PRI population was not included due to the single locality from Russia; (ii) Longitude, populations were divided into two groups (West and East) by the central line of the Tsushima Strait (130°E); (iii) Biogeographic basins (Wang, 1999), populations were divided into four groups: Japan–Pacific coast (JPC), Sea of Japan (SOJ), East China Sea including Yellow–Bohai Sea (ECS), and South China Sea (SCS); (iv) Marine provinces (Longhurst, 2007), populations were divided into two groups: the Kuroshio Current province (KURO) and China Sea Coastal province (CHIN). The significance (p = 0.05) of the fixation indices was calculated based on 10⁴ permutations. Pairwise genetic differentiations (F_st) and Nei’s genetic distance among lineages were calculated using Arlequin and POPGENE v.1.32 (Yeh et al., 1999). In addition, pairwise genetic differentiations (F_st) and gene flow (N_m) among populations were also calculated in Arlequin. Based on the estimated substitution rate of mtDNA cox1 in the class Florideophyceae (Le Gall and Saunders, 2010; Bringloe and Saunders, 2019), we set a range of 0.7–1.0%/Ma to broadly calculate divergence times between cox1 lineages in A. vermiculophyllum in the native northwest Pacific.

S-DIVA (Statistical Dispersal-Vicariance Analysis) implemented in RASP v.3.1 (Yu et al., 2015) was used to infer the phylogeographic history of A. vermiculophyllum through statistically evaluating the alternative ancestral ranges at each node in a phylogenetic tree. The tree data set and condensed tree produced by BEAST and TreeAnnotator were loaded in RASP with the outgroup removed and then saved as “trees file.” Four main biogeographic regions (Wang, 1999) were defined based on geographic distribution of cox1 haplotypes: (A) Japan–Pacific coast, (B) Sea of Japan, (C) East China Sea including Yellow–Bohai Sea, and (D) South China Sea. The “Distribution file” was produced based on the distribution of cox1 haplotypes through all defined biogeographic regions. “Trees file” and “distribution file” were loaded in RASP to perform the S-DIVA analysis. The topology and Most Likely States (MLS) functions were used to infer the probability of ancestor distribution at each node of the phylogeny.

The length of the cox1 alignment was 1,214 nucleotides (nt) and consisted of 611 sequences. A total of 45 polymorphic sites (3.71%) were detected (Supplementary Table S2), including 28 singletons (2.31%) and 17 parsimony informative sites (1.40%). Thirty-three haplotypes were defined (GenBank accession numbers: MN817547-MN817579), of which 26 were restricted to a single population and seven were found in two or more populations (Supplementary Table S1). The two dominant haplotypes H1 (90 specimens, 15%) and H3 (339 specimens, 55%) were found in 11 (23%) and 28 (58%) of the populations, respectively (Supplementary Table S1). There were 28 populations harboring a single haplotype and six populations had more than two haplotypes (maximum four) (Supplementary Table S1). Haplotype diversity (total Hd = 0.6632) ranged from 0.0000 to 0.7905, and nucleotide diversity (total π = 0.2473 × 10^–2) varied from 0.0000 to 0.3766 × 10^–2 (Supplementary Table S1). The east of Jeju Island population (JEE, 31 in Figure 1), Korea had the highest haplotype diversity (Hd = 0.7905), and the south of Jeju Island population (JES, 32 in Figure 1), Korea had the highest nucleotide diversity (π = 0.3766 × 10^–2) (Supplementary Table S1). At a regional scale, populations from Japan showed much higher genetic diversity (133 specimens, Hd = 0.7236, π = 0.3275 × 10^–2) than those from Korea (213 specimens, Hd = 0.5282, π = 0.1361 × 10^–2), China (262 specimens, Hd = 0.5265, π = 0.2100 × 10^–2) and Russia (three specimens, Hd = 0.0000, π = 0.0000) (Supplementary Table S1). The value of neutrality test of the A. vermiculophyllum populations in the northwest Pacific is negative (Tajima’s D = −1.4414, p < 0.05; Fu’s Fs = -10.1928, p < 0.05) (Supplementary Table S1).

Phylogenetic analyses of the 33 A. vermiculophyllum cox1 haplotypes with two outgroup taxa resulted in five strongly supported A. vermiculophyllum clades (lineages 1–5 in Figure 2) in both the BI and ML. BI and ML trees presented highly concordant phylogenetic structure. Haplotype 10 differed from haplotype 32/33 by eight bases, so haplotype 10 was regarded as an independent genetic lineage (Figure 3). Another four diverged lineages were supported with high values for main nodes (BI: posterior probabilities = 0.71∼0.99, ML: bootstrap values = 78∼96) except for lineage 5 in the ML (Figure 2). The haplotype network (Figure 3) revealed the same five distinct lineages in A. vermiculophyllum which were basically in accordance with phylogenetic trees. Both lineages 2 and 5 showed significant star-like structure, and the neutrality test of them was negative (Figure 3 and Table 1). In addition, there were 3–11 substitutions detected among the five lineages (e.g., lineage 2 differed by three bases from lineage 5, and lineage 1 differed by 11 bases from lineage 4) (Figure 3). Positions 227, 329, 377, 398, 443, 485, 722, 767, 860, 945, 989, 1040, 1106, and 1154 of the cox1 alignment showed nucleotide substitutions unique to certain lineages (Supplementary Table S2). For example, at the 398^th site, nucleotide A corresponded to lineage 5, while nucleotide G corresponded to the other four lineages. Similar distribution patterns of nucleotides were found at other sites. In addition, the optimal K value (K = 5) was detected by genetic clustering analysis (Supplementary Figure S1), which was consistent with phylogenetic results and biogeographical distribution of haplotypes (Figures 1, 2).

Mapping the distribution of the five lineages showed a basic biogeographic pattern (Figure 1): Lineage 1 consisting of five haplotypes was mainly distributed in the south of Japan–Pacific coast and the Sea of Japan. Lineage 2 consisting of six haplotypes was dominated by H1 and mainly distributed in the north of Japan–Pacific coast and the Sea of Japan. Lineage 3 was only distributed in one population in the south of Japan–Pacific. Lineage 4 consisting of two haplotypes was distributed in three southern Chinese populations. The most widely distributed Lineage 5 consisting of 19 haplotypes was dominated by H3 and mainly distributed in the East China Sea, the Yellow–Bohai Sea, and the Sea of Japan. Most populations (42/48) only harbored one lineage (Figure 1). Lineage 1 exhibited higher levels of genetic diversity (Hd = 0.7988, π = 0.2140 × 10^–2) than the others (Table 1).

AMOVA revealed that most genetic variations occurred among populations within groups (67–75%) except the grouping criterion of biogeographic basins (61% of the variation occurred among groups), and all measures were statistically significant (p < 0.001) (Table 2). Pairwise F_st estimates exhibited high levels of genetic differentiation among lineages (F_st > 0.8000, p < 0.001), and Nei’s genetic distance varied between 0.0726 and 0.2716 (Table 3). Population-level F_st values also revealed high levels of genetic differentiation and limited connection among populations (70.8% of F_st > 0.25, p < 0.05; 11.3% of N_m > 1) (Supplementary Table S3). In particular, the geographically proximate populations (13, 14, and 15) in the south of Japan–Pacific were highly differentiated (F_st = 1.00, p < 0.001; N_m = 0.0000). The populations 11, 12, 14, and 15 from Japan, the populations 22 and 29 from Korea, and the populations SS and ZP from China were highly differentiated from the others (F_st = 0.4209–1.0000, p < 0.05) (Supplementary Table S3). Sequence divergence among five lineages varied from 0.41 to 1.03%, and the estimated divergence time between them dating to 0.21–0.74 Ma.

S-DIVA analysis of mtDNA cox1 haplotypes suggested that current A. vermiculophyllum populations originated from a wide ancestral region encompassing the Japan–Pacific coast, East and South China Sea (node 65, ACD in Figure 4), associated with multiple dispersal and vicariance events (Figure 4). Lineage 1 originated in the Japan–Pacific coast and the Sea of Japan (node37, AB); lineage 2 originated along the Japan–Pacific coast (node 42, A); Lineages 3 and 4 were isolated due to vicariance events from the common ancestor distributed along Japan–Pacific coast and the South China Sea (node 44, AD); lineage 5 originated from the East China Sea (node 63, C). A total of 12 dispersal events and eight vicariance events were detected (Figure 4).

In this study, the detected overall haplotype (Hd = 0.6632) and nucleotide diversity (π = 0.3766 × 10^–2) are similar to previous results reported in A. vermiculophyllum in the northwest Pacific by Kim et al. (2010) (Hd = 0.731, π = 0.329 × 10^–2) and Liu et al. (2016) (Hd = 0.589, π = 0.318 × 10^–2), respectively. Such high haplotype diversity and low nucleotide diversity patterns suggest that A. vermiculophyllum in the native northwest Pacific may have experienced rapid population growth over a short period because sudden population expansion cannot support this species to have sufficient time to accumulate enough nucleotide mutations (Grant and Bowen, 1998). The haplotype network of A. vermiuclophyllum (e.g., lineages 2 and 5 in Figure 3) hints to recent population expansion as well because new haplotypes derived from ancestral haplotypes (H3 and H1) usually formed a star-like topology during population expansion (Posada and Crandall, 2001). This conjecture of population expansion is also concordant with the neutrality test (Tajima’s D = −1.4414, p < 0.05; Fu’s Fs = −10.1928, p < 0.05; Supplementary Table S1).

The high haplotype diversity of A. vermivulophyllum in the northwest Pacific may be ascribed to high environmental heterogeneity (Nei, 1987). Discontinuous habitat (e.g., Fucus vesiculosus, Zardi et al., 2013) and variation of salinity (e.g., Cladophora, Hayakawa et al., 2012) have been reported to impact genetic diversity and differentiation of intertidal seaweeds. A. vermiculophyllum is distributed in both estuaries and marine coast in the northwest Pacific. Regional habitat discontinuity and heterogeneity (e.g., sandy beach, rocky reef, gravel and mud beach, Supplementary Table S1) may have caused population-level genetic differentiation. In addition, A. vermiculophyllum is well adapted to salinity variation of 5–60‰ (Yokoya et al., 1999). Such salinity gradients may have imposed strong selection on allele frequencies (e.g., Furcellaria lumbricalis, Kostamo et al., 2012) and apomictic reproduction and phenotypic variation (e.g., Ulva prolifera, Ogawa et al., 2014), ultimately leading to genotypic differentiation.

We detected five genetic lineages in A. vermiculophyllum in the native northwest Pacific. Three of them (lineages 1, 2, and 5) showed similar distribution ranges as reported by Krueger-Hadfield et al. (2017). However, the ∼35°N segregation based on nucleotides at sites 945 and 1040 is not as clear as in Krueger-Hadfield et al. (2017). Although haplotypes with a T at these positions all occur north of ∼35°N (e.g., lineage 2), haplotypes with a C were found on both sides of this latitudinal divide. Twelve additional unique nucleotide substitutions also correspond to lineage differentiation, respectively (Supplementary Table S2). Our data suggest that a gradual nucleotide change in lineage distribution happened at 35°N, and that nucleotide substitutions at a few sites are not sufficient to delineate genetic lineages and their geographic extent as described by Krueger-Hadfield et al. (2017).

Sea levels in the northwest Pacific dropped by 120–140 m during the glacial maxima in the late Quaternary (e.g., the Last Glacial Maximum, Lambeck et al., 2002), leading to the isolation of the Sea of Japan and a shift of the coastline of the South China Sea and the East China Sea (and the presumed disappearance of the Yellow–Bohai Sea) (Wang, 1999). This coastline configuration imposed significant impacts on distribution range and genetic diversity of marine species (Cheang et al., 2010; Ni et al., 2014; Hu et al., 2015). Our S-DIVA analysis illustrated that the current distribution range of A. vermiculophyllum stem from the ancestral relics that persisted historical glaciation periods and experienced multiple dispersal events (Figure 4). It should be noted that the sister relationships between the five lineages have not been resolved phylogenetically. Sequencing multiple genes or complete mitochondria of specimens from each lineage may help to address this issue. Phylogeographic studies have proposed Japan–Pacific coast (Hu et al., 2011, 2015, 2017), Okinawa trough (Wang, 1999; Hu et al., 2015), and northern South China Sea (Hu et al., 2018a) as glacial refugia in the northwest Pacific. The locations of these glacial refugia broadly coincide with the distribution of five lineages in A. vermiculophyllum in the native northwest Pacific. For example, lineage 2 originated from the north of Japan–Pacific coast and spread southward into the Sea of Japan driven by the Oyashio Current and the Tsushima Warm Current. Similarly, lineages 4 and 5 originated from the South China Sea and the Okinawa trough, respectively, and colonized the South China Sea coast (lineage 4) and the Yellow–Bohai Sea (lineage 5) driven by the South China Sea Warm Current, the Yellow Sea Warm Current and the China Coastal Current (Figures 1, 4). Mitochondrial lineages 2/3/4, despite their geographic distance in the northwest Pacific, exhibit close genetic relationships, which likely resulted from the persistence of the most recent common ancestor in isolated refugia during the ice age and subsequently multiple dispersal and vicariance events (node 45 in Figure 4). Therefore, ancestral persistence in multiple scattered glacial refugia during the late Quaternary ice ages is likely the major phylogeographic process contributing to extant lineage differentiation and distribution.

Geographically, we found two distinct genetic breaks (Figure 1). The first break occurred in the south-central coast of China, which exactly corresponds to the boundary between the distribution range of lineages 4 and 5 (Figure 1). As illustrated by Hu et al. (2018b), the genetic break may have resulted from the physical barrier of the Yangtze river diluted water. The diluted water is characterized by profoundly decreased seawater salinity and pH, which can affect the trans-regional dispersal of intertidal marine organisms (Wang et al., 2003) and hence population-level genetic connectivity. It has been reported that the growth of A. vermiculophyllum is inhibited by salinity fluctuation, and the greater the salinity fluctuation, the more obvious the inhibition (Nejrupa and Pedersena, 2012). In northern China, spore release of A. vermiculophyllum occurs mainly in summer, especially from June to July (Li and Li, 1986). The Yangtze river diluted water, which is influenced by southwestward monsoon in summer, deflected to the northeast paralleling with the Taiwan Warm Current, consequently caused the deflection of southward China Coastal Current and prevented the northern populations spreading to the south (Dong et al., 2012). Alternatively, the genetic break may have occurred near the Taiwan Strait as the reported genetic split between the brown algae Colpomenia peregrina and Colpomenia sinuosa (Song et al., 2019). Such a break may be the result of the isolation of the Taiwan Strait land bridge caused by the sea level drop during glaciation periods (Ota, 1998), leading to the isolation of species/populations on both sides of the land bridge which have also long been defined as marine ecoregions 52 and 113 (Spalding et al., 2007). To test this hypothesis, we need to sample more populations from the region between populations SS (45 in Figure 1) and ZP (46 in Figure 1). It should be noted that lineage 4 has close relatives (lineages 2 and 3) that exist further north than lineage 5. Therefore, gene flow probably had occurred in the past and the Yangtze river was not the primary vicariance event that separated lineages 4 and 5.

Another genetic break was detected between lineages 1 and 2 in the central Japan–Pacific coast (Figure 1) which corresponds to the boundary of the marine ecoregions 48 and 51 (Spalding et al., 2007), suggesting a long-term genetic isolation despite their geographic proximity. This differentiation pattern is similar to results found in Sargassum horneri (Turner) C. Agardh and Sargassum fusiforme (Harv.) Setchel in the northwest Pacific (Hu et al., 2011, 2017). The central Japan–Pacific coast is characterized by the convergence of the cold Oyashio Current and the warm Kuroshio Current (Figure 1). The interplaying current systems, especially the coastal upwelling, can create distinct physical barrier for genetic exchange, particularly for the passive drifting spores and fragments of A. vermiculophyllum (Yu and Kim, 2018). However, the current and upwelling regime may not have served as the primary vicariance separating lineages 1 and 2 because they are only distantly related. Lineage 2 is specifically restricted to the north of the Japanese Archipelago and east coast Korean Peninsula, and the underlying phylogeographic processes probably correlates with contemporary ocean currents, but more evidence is needed with wide-range sampling.

Agarophyton vermiculophyllum is an important economic and ecological species that deserves to be protected in the situation of increasing biodiversity destruction in the northwest Pacific (Armisen, 1995; Abreu et al., 2011; Tanaka et al., 2012). Our mtDNA-based phylogeographical study provides a basis for how to conserve A. vermiculophyllum populations in the wild. We detected the highest population genetic diversity and number of genetic lineages in the Japanese Archipelago (Figure 1 and Supplementary Table S1). In particular, this area is rich in unique and endemic haplotypes, and the F_st values between adjacent populations are much higher than those between the distant Korean and Japanese populations in the south of Japan–Pacific coast, reflecting a potential habitat fragmentation (Rossiter et al., 2007; Hu et al., 2017). In addition, gene flow between populations in this area is extremely low (N_m = 0.0000), indicating very low population connectivity. In order to prevent the loss of local genetic diversity caused by a genetic drift, the Japanese Archipelago can be considered as a special conservation priority. In addition, the five genetic lineages with significant genetic differentiation were isolated since the late Quaternary ice ages with significant genetic differentiation. Particularly, lineage 5 has the highest genetic diversity and is reciprocally monophyletic with the other four lineages; it thus should be taken as an independent Evolutionary significant unit (ESU) to be preserved to maintain commercial access and ecological sustainability (Moritz, 1994; Hu et al., 2017).