MtDNA-Based Phylogeography of the Red Alga Agarophyton vermiculophyllum (Gigartinales, Rhodophyta) in the Native Northwest Pacific

The repeated transgression and regression of coastlines mediated by the late Quaternary glacial–interglacial cycles make the northwest Pacific a hot spot to study marine speciation and population diversity. The red alga Agarophyton vermiculophyllum is an ecologically important species native to the northwest Pacific, capturing considerable research interest due to its wide-range invasiveness in Europe and North America. However, the knowledge of phylogeographic structure and intraspecific genetic diversity across the entire native range was still scarce. Here, we used 1,214-bp of mitochondrial cytochrome c oxidase subunit 1 (cox1) to explore phylogeographic patterns, lineage structure, and population genetic differentiation of 48 A. vermiculophyllum populations in the northwest Pacific. Our DNA data revealed overall high haplotype diversity and low nucleotide diversity and five phylogeographically structured genetic lineages that diverged significantly from each other. S-DIVA analysis showed the ancestors of A. vermiculophyllum originating from multiple areas encompassing the Japan–Pacific coast, East and South China Seas. This combined evidence indicates that A. vermiculophyllum might have survived in multiple scattered glacial refugia during the late Quaternary climate oscillations in the northwest Pacific. Such knowledge may help to better understand how palaeoclimate interacted with contemporary environments to contribute to intraspecific genetic variation and provide a new perspective for conserving natural resource of A. vermiculophyllum in the northwest Pacific.

section_content: The northwest Pacific spans three major marginal seas, the Sea of Japan, the East China Sea (including Yellow-Bohai Sea), and the South China Sea (Figure 1 ). In the late Quaternary, global sea level fell 120-140 m during glacial maxima and thus led to the exposure of the marginal seas (Lambeck et al., 2002). Consequently, the Yellow-Bohai Sea was completely exposed, the East China Sea shrank eastward to the Okinawa trough, the South China Sea became a semi-closed sea, and the Sea of Japan with a shallower bedrock became an almost enclosed sea (Figure 1 ; Wang, 1999; Oba and Irino, 2012). During the interglacial periods, the three marginal seas were inundated by rising sea levels. These repeated transgression-regression cycles, together with long-distance dispersal mediated by ocean currents and recent anthropogenic activities, imposed great impacts on phylogeographic patterns and population genetic differentiation in this area (Cheang et al., 2010; Wang et al., 2016). The northwest Pacific has therefore been proposed as an ideal model for studying how sea level oscillations, caused by repeated glaciations, structured distribution range and spatial population genetic diversity and marine phylogeography (Ni et al., 2014). In the northwest Pacific, the reduction of marginal seas during the late Quaternary is thought to have led to a large decrease of marine species, while some ancestral relict populations may have survived in isolated refugia. The discontinuity of these refugia made an essential contribution to population-level isolation and differentiation (Ni et al., 2014). Phylogeographic studies of a variety of marine organisms have proposed several marine glacial refugia in the northwest Pacific, which scattered from the South China Sea (Hu et al., 2018a) to Taiwan island (Lee et al., 2012), the Okinawa trough (Hu et al., 2015), the Sea of Japan (Wang, 1999), Jeju island (Lee et al., 2012), and the Japanese Pacific coast (Hu et al., 2011 (Hu et al.,, 2015 (Hu et al.,, 2017)). When sea levels rose after the Last Glacial Maximum (LGM), the surviving populations in marine refugia expanded to colonize new habitats, eliminating population differentiation caused by isolation in different refugia (Xue et al., 2014). Ocean currents form another important driving force of post-glacial population colonization. The dominant ocean current systems in the northwest Pacific constitute the Kuroshio Current and its branches (the Yellow Sea Warm Current and the Tsushima Warm Current), the China Coastal Current, the South China Sea Warm Current, and the Oyashio Current (Figure 1 ). In the Yellow-Bohai Sea, for example, the Yellow Sea Warm Current originating from the Okinawa trough has been reported to drive genetic homogeneity in some marine organisms (Hu et al., 2015; Li et al., 2017; Liu et al., 2018). 

The red alga Agarophyton vermiculophyllum (Ohmi) Gurgel, J. N. Norris et Fredericq, is native to the northwest Pacific (Gurgel et al., 2018). It was originally described as Gracilariopsis vermiculophylla Ohmi from the Akkeshi Bay, Hokkaido, Japan, but was transferred a decade later to the genus Gracilaria Greville and only recently was transferred to a new genus Agarophyton Gurgel, J. N. Norris et Fredericq, together with two other former Gracilaria species (Gurgel et al., 2018). Its life cycle consists of an alternation of isomorphic haploid and diploid stages. It is considered an important economic and ecological species because it is one of the main raw materials for extracting agar (Armisen, 1995) and used as a biofilter in seawater eutrophication restoration (Abreu et al., 2011). 

The past few decades witnessed an invasion of A. vermiculophyllum in Europe and North America (Krueger-Hadfield et al., 2017). The invasion of A. vermiculophyllum has captured widespread research interest, particularly its impacts on coastal ecosystem and adaptation to environmental variables (Davoult et al., 2017; Sotka et al., 2018). Some studies indicate that strong asexual reproduction, high tolerance to abiotic stresses, low palatability, and high resistance to epiphytes account for the rapid invasion, establishment, and expansion in new habitats (Hu and Lopez-Bautista, 2014). Kim et al. (2010) and Krueger-Hadfield et al. (2017) speculated the Sea of Japan and northeastern Japan as the donor areas of invasive A. vermiculophyllum populations in Europe and North America, possibly mediated by the exported oyster Crassostrea gigas Thunberg, 1793 from Japan, but these studies lacked adequate sampling from its native range. 

Molecular studies of A. vermiculophyllum in its native range are still scarce. Kim et al. (2010) and Gulbransen et al. (2012) compared genetic diversity based on mitochondrial cox1 (cytochrome c oxidase subunit 1) and found that A. vermiculophyllum in the native northwest Pacific had substantially higher genetic diversity and haplotype richness compared to the invasive populations in Europe and North America. Krueger-Hadfield et al. (2017) detected three mitochondrial lineages of A. vermiculophyllum in the northwest Pacific and high genetic differentiation among proximal geographic populations using microsatellite genotyping. They found two cox1 genetic groups that differed by two C/T substitutions which showed a segregated distribution pattern north and south of ∼35 • N. Recently, Hu et al. (2018b) reported a significant genetic break in A. vermiculophyllum along the coast of China using Amplified Fragment Length Polymorphism data. Nevertheless, they did not yet thoroughly understand phylogeographic structure and potential cryptic lineage diversity in A. vermiculophyllum across its native rang in the northwest Pacific. 

In the northwest Pacific, A. vermiculophyllum is distributed from the far east of Russia in the north to tropical Guangxi, China, and Vietnam in the south (Xia, 1999; Skriptsova and Choi, 2009). In this study, we expanded the sampling area of A. vermiculophyllum in its native range in the northwest Pacific compared to previous studies to (i) quantify phylogeographic structure and genetic variation of A. vermiculophyllum in the northwest Pacific and (ii) detect the diverged genetic lineages and reconstruct their putative phylogeographic histories. We are especially interested in testing the hypothesis proposed by recent phylogeographical studies (Hu et al., 2011; Dong et al., 2012; Ni et al., 2014; Li et al., 2017) in which sea level drops during the Quaternary ice ages and oceanic currents interactively contributed to geographical distribution and genetic differentiation of coastal marine organisms in the northwest Pacific. This knowledge may enhance our understanding of S1. Colors corresponds to the lineages indicated in the phylogenetic tree (Figure 2 ) and haplotype network (Figure 3 ). The black areas represent the elevation >0 m, the gray areas represent the depth of water <130 m, and the light blue areas represent the depth of water >130 m (sea levels during the last glacial maximum). KC, Kuroshio Current; TWC1, Tsushima Warm Current; YSWC, Yellow Sea Warm Current; CCC, China Coastal Current; TWC2, Taiwan Warm Current; SCSWC, South China Sea Warm Current; OC, Oyashio Current. 

the distribution shifts and genetic exchange patterns of marine species in this area and provide action guidelines for conserving natural resources of A. vermiculophyllum.
section_name: INTRODUCTION

section_name: MATERIALS AND METHODS

section_content: 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.
section_name: Sample Collection

section_content: 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).
section_name: DNA Extraction, PCR Amplification and Sequencing

section_content: 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 7 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) 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 4 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 4 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.
section_name: Molecular Diversity, Phylogeny, and Genetic Structure

section_content: 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.
section_name: Historical Phylogeography

section_name: RESULTS

section_content: 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 ).
section_name: Genetic Diversity

section_content: 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 Positions 227, 329, 377, 398, 443, 485, 722, 767, 860, 945, 989,, 1106 Positions 227, 329, 377, 398, 443, 485, 722, 767, 860, 945, 989,, 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.
section_name: Phylogenetic and Phylogeographic Structure

section_content: 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 : 2, [5] [6] [7] [8] [9] 11, 12, 14, 15, [30] [31] [32] SOJ: 1, 3, 4, 10, 13, 16, 17, 22; [23] [24] [25] [26] [27] [28] [29] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] 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 ). (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 populationlevel 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).
section_name: Historical Phylogeography

section_content: 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 (Hu et al.,, 2015 (Hu et al.,, 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 (Hu et al.,, 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.
section_name: Phylogeographic History of A. vermiculophyllum in the Northwest Pacific

section_content: 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).
section_name: Conservation Insights for A. vermiculophyllum

section_name: ACKNOWLEDGMENTS

section_content: We thank Dr. Zhong-Min Sun and Mr. Ruo-Yu Liu for providing assistance in field collections.
section_name: ACKNOWLEDGMENTS

section_name: FUNDING

section_content: This study was supported by the National Natural Science Foundation of China ( 41761144057 and 31971395 ), the Sino-German Centre for Research Promotion ( GZ1375 ), the Strategic Priority Research Program of the Chinese Academy of Sciences ( XDB42000000 ), and the Thailand Research Fund ( RDG6130002 ).
section_name: FUNDING

section_name: DATA AVAILABILITY STATEMENT

section_content: The datasets generated for this study can be found in NCBI GenBank, NCBI Accession No. 's MN817547-MN817579.
section_name: DATA AVAILABILITY STATEMENT

section_name: AUTHOR CONTRIBUTIONS

section_content: Z-MH conceived the project. K-LZ conducted the molecular experiments, data collection, analysis, and writing. X-HS assisted in data analysis and writing. H-GC, SS, and FW provided the partial samples and assisted data interpretation. SD, D-LD, and Z-MH commented the results and modified the manuscript to final version.
section_name: AUTHOR CONTRIBUTIONS

section_content: The Supplementary Material for this article can be found online at: https://www. frontiersin. org/articles/10. 3389/fmars. 2020. 00366/full#supplementary-material
section_name: SUPPLEMENTARY MATERIAL

section_content: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
section_name: Conflict of Interest: