Historical Connectivity and Demography of the Ferocious Reef Crab, Eriphia ferox (Crustacea; Eriphiidae), Demonstrate That Taoyuan Algal Reef Is an Essential Population Source Along the East Taiwan Strait

The east Taiwan Strait is largely fringed by sandy and muddy habitats. However, a massive algal reef made of crustose coralline algae has been found along the coast off Taoyuan city in northwestern Taiwan. The porous structure of Taoyuan Algal Reef harbors high abundance and diversity in marine organisms, including the ferocious reef crab, Eriphia ferox. Such a pivotal geographic location and unique ecological features make Taoyuan Algal Reef a potential stepping stone connecting biotic reefs in the east Taiwan Strait, South China Sea to the south, and even the high latitude of Japan to the north. In this study, we examined the population connectivity and historical demography of E. ferox by analyzing mitochondrial cytochrome oxidase I (COI) fragments of 317 individuals sampled from 21 localities in the northwestern Pacific. Our analyses of haplotype network and pairwise FST comparisons revealed a lack of phylogeographical structure among E. ferox populations, implying the existence of a migration corridor connecting the South and East China Seas through the east Taiwan Strait. Multiple lines of evidence, including significant values in neutrality tests, unimodally shaped mismatch distributions, and Bayesian skyline plots elucidated the rapid population growth of E. ferox following the sea-level rise after Last Glacial Maximum (ca. 2–10 Ka). Such demographic expansion in E. ferox coincided with the time when Taoyuan Algal Reef started to build up around 7,500 years ago. Coalescent migration analyses further indicated that the large and continuous E. ferox population exclusively found in Datan Algal Reef, the heart of Taoyuan Algal Reef, was a source population exporting migrants both northward and southward to the adjacent populations. The bidirectional gene flow should be attributed to larval dispersal by ocean currents and secondary contact due to historical population expansion. Instead of serving as a stepping stone, our results support that Taoyuan Algal Reef is an essential population source for biotic reef-associated species along the east Taiwan Strait, and highlight the importance of conserving such a unique ecosystem currently threatened by anthropogenic development.

The east Taiwan Strait is largely fringed by sandy and muddy habitats. However, a massive algal reef made of crustose coralline algae has been found along the coast off Taoyuan city in northwestern Taiwan. The porous structure of Taoyuan Algal Reef harbors high abundance and diversity in marine organisms, including the ferocious reef crab, Eriphia ferox. Such a pivotal geographic location and unique ecological features make Taoyuan Algal Reef a potential stepping stone connecting biotic reefs in the east Taiwan Strait, South China Sea to the south, and even the high latitude of Japan to the north. In this study, we examined the population connectivity and historical demography of E. ferox by analyzing mitochondrial cytochrome oxidase I (COI) fragments of 317 individuals sampled from 21 localities in the northwestern Pacific. Our analyses of haplotype network and pairwise F ST comparisons revealed a lack of phylogeographical structure among E. ferox populations, implying the existence of a migration corridor connecting the South and East China Seas through the east Taiwan Strait. Multiple lines of evidence, including significant values in neutrality tests, unimodally shaped mismatch distributions, and Bayesian skyline plots elucidated the rapid population growth of E. ferox following the sea-level rise after Last Glacial Maximum (ca. 2-10 Ka). Such demographic expansion in E. ferox coincided with the time when Taoyuan Algal Reef started to build up around 7,500 years ago. Coalescent migration analyses further indicated that the large and continuous E. ferox population exclusively found in Datan Algal Reef, the heart of Taoyuan Algal Reef, was a source population exporting migrants both northward and southward to the adjacent

INTRODUCTION
The Taiwan Strait, bounded between southeastern China and Taiwan, is a shallow passage bridging the South China Sea (SCS) and the East China Sea (ECS). Situated at the tectonic and geographic junction along the coast of eastern Eurasia, historical closure and opening of the Taiwan Strait were affected profoundly by sea-level fluctuations during the last glacial cycles. The interplay of paleoclimate and distinct topographical features of the Taiwan Strait has subsequently shaped the distribution patterns and genetic diversity of extant marine organisms in the northwestern Pacific (Ni et al., 2014 and references therein).
During glacial periods throughout the Pleistocene (2.6 Ma-11.7 ka), the recurrent desiccation and exposure of land mass in the Taiwan Strait culminated in transient disconnections between SCS and ECS basins (Wang, 1999;Kimura, 2000;Ludt and Rocha, 2015). Consequently, gene flow between the two marginal seas was blocked and divergence in marine organisms occurred. Distinct lineages were found to be concordant with the two sea basins, which should have served as independent refugia during glaciations (Ni et al., 2014;Wang et al., 2017;Qu et al., 2021 and references therein). As glaciers receded at the end of the Last Glacial Maximum (LGM,(17)(18)(19), sea level rose drastically and flooded the entire Taiwan Strait, establishing the coastlines and ocean current regimes that we know today . The emerging waterway conceptually restores oceanic connectivity, and should facilitate genetic homogenization in marine organisms along the Taiwan Strait. Surprisingly, studies of moon turban snail Lunella granulata (Gmelin, 1791) and Portuguese oyster Magallana (formerly Crassostrea) angulata (Lamarck, 1819) revealed a significant subdivision in populations between the Chinese and Taiwanese coasts of the Taiwan Strait (Chiu et al., 2013;Hsiao et al., 2016). As our current understanding of marine phylogeography in the northwestern Pacific exclusively relies on studies conducted along the Chinese coast (reviewed in Ni et al., 2014), such cross-strait divergence highlights the necessity to explore the role of the uninvestigated west coast of Taiwan.
The west coast of Taiwan is mainly composed of sandy/muddy habitats (Chen and Shashank, 2009). Hard substrates constructed by reef-building corals appear only at the northern/southern tips and east coast of Taiwan (Dai and Cheng, 2020). The only biotic reef along the west coast of Taiwan is a massive algal reef built by crustose coralline algae (CCA) in the intertidal zone to a depth of around 10 m along a 27 km stretch of the coastline off Taoyuan city (Liou et al., 2017). The porous algal reef has an ancient origin, dated ca. 7,500 BP (Dai et al., 2009). Its complex structures provide shelter for diverse marine organisms, including crustaceans, polychaetes, sipunculans, and the largest population of the critically-endangered caryophylliid coral Polycyathus chaishanensis (Lin et al., 2012) (Liou et al., 2017;Kuo et al., 2020) in Taiwan. In addition, Taoyuan Algal Reef is reported to share specific marine species like Oulastrea crispata (Lamarck, 1816) and Epinephelus quoyanus (Valenciennes, 1830) with Penghu archipelago, which comprises as many as 90 islands and islets in the southeastern Taiwan Strait. These organisms are not found elsewhere in the Taiwan Strait (Chen, 2017). Its ancient origin and pivotal geographic location, supplemented with the relatively high diversity of marine benthos and species connectivity with the Penghu Islands, make Taoyuan Algal Reef a potential stepping stone connecting biotic reefs in the east Taiwan Strait (Chen, 2017;Kuo et al., 2020).
To examine the hypothesis of Taoyuan Algal Reef as a stepping stone, we investigated the genetic structure of the ferocious reef crab Eriphia ferox (Koh and Ng, 2008), a species recognized as Eriphia smithii (MacLeay, 1836) until 2008 (Koh and Ng, 2008). The relatively newly described E. ferox has a widespread distribution along the hard substrates of the northwestern Pacific, ranging from southern Japan in the north to the southern SCS in the south (Koh and Ng, 2008). In 2019, a significantly large and continuous population of E. ferox was found at Datan , which is located within the Taoyuan Algal Reef. These biogeographical and ecological characteristics make E. ferox an ideal model organism to test the stepping stone hypothesis. We used the mitochondrial (mt) cytochrome oxidase I (COI) fragment to elucidate the genetic structure, direction of gene flow, and historical demography and connectivity of E. ferox. The findings of the present study should broaden our current understanding of marine phylogeography in the northwestern Pacific.

Sampling and DNA Sequencing
In the field, crabs were captured for tissue sampling and released afterward. From each crab, one pereiopod was removed and preserved in 95% ethanol for DNA extraction later. A total of 317 samples were collected from southern Japan, northern/northwestern Taiwan (mainland), the Penghu archipelago (Taiwan), Malaysia, and Singapore (Table 1). Genomic DNA was extracted from pereiopod muscle tissue using a Tissue and Cell Genomic DNA Purification Kit (GeneMark, Taiwan) according to the company's protocol. A fragment containing the COI gene was amplified with the E. ferox specific primer set ErifeCOI_F (forward primer, 5 -GTA GGA ACT TCA CTA AGA CTA ATT ATT CG-3 ) and ErifeCOI_R (reverse primer, 5 -GCA GGG TCA AAG AAG GAT G-3 ), which was manually designed in this study. The COI sequence (HM638034.1) by Lai et al. (2014) was used as reference. The volume of the PCR mixture was 25 µL, containing 0.5-1.0 µL of genomic DNA template, 0.5 µL of each primer (10 µM), 12.5 µL of Taq DNA Polymerase Master Mix RED (Ampliqon, Denmark), and 10.5 µL of ddH 2 O. PCR was performed under the following conditions: initial denaturation at 95 • C for 3 min, followed by 40 cycles at 94 • C for 30 s, 50-55 • C for 40 s, 72 • C for 30 s, and a final extension at 72 • C for 5 min. The resulting PCR products were sent to the company Genomics (Taipei, Taiwan) for Sanger sequencing.

Data Analyses
Raw sequences were inspected, trimmed and assembled with the software Geneious version R10 (Biomatters, New Zealand). Cleaned consensus sequence of each sample was imported to MEGA X (Kumar et al., 2018) for sequence alignment using the built-in MUSCLE algorithm. Genetic diversity was evaluated for the number of haplotypes (n h ), haplotype diversity (h), and nucleotide diversity (π) for each sampling site using DnaSP 6.12 (Rozas et al., 2017). The haplotype network was illustrated with the median-joining algorithm in POPART (Leigh and Bryant, 2015). Pairwise F ST comparisons were calculated between sampling sites, with statistical significance assessed by performing 10,000 permutations using Arlequin 3.5.2 (Excoffier and Lischer, 2010). Statistical significance of F ST estimates was subsequently adjusted by Benjamini and Hochberg's FDR correction (Benjamini and Hochberg, 1995) in R version 3.6.1 (R Core Team, 2013). To investigate regional demographic history, a total of five groups were defined for the analyses: Japan, Taiwan (pooled northern/northwestern Taiwanese populations), Penghu (pooled Penghu populations), Malaysia (pooled Malaysian populations), and Singapore. Signatures of historical demography were inferred from neutrality tests, mismatch distributions, and Bayesian Skyline Plots. For neutrality tests, Tajima's D and Fu's Fs were calculated by Arlequin with 10,000 permutations (Tajima, 1989;Fu, 1997), while Ramos-Onsins and Roza's R 2 index (Ramos-Onsins and Rozas, 2002) were calculated using DnaSP. For every geographic population, mismatch distribution was established to test against a model of exponential population growth (Rogers and Harpending, 1992). The validity of the sudden demographic expansion model was examined in Arlequin by computing the Harpending's raggedness index (rg, Harpending, 1994) and sum of square deviations (SSD) between the observed and expected mismatch distributions. The dynamics of effective population size (N e ) through time were visualized with Bayesian Skyline Plots (BSPs), which were performed using BEAST 2.6.5 (Bouckaert et al., 2019). The HKY model was selected as the best fit model by jModelTest 2.1.10 (Darriba et al., 2012), and a mutation rate of 4.87% TABLE 1 | Sampling sites, sample size (n), number of haplotypes (n h ), haplotype diversity (h ± SD), and nucleotide diversity (π ± SD) from samples of Eriphia ferox collected along the west coast of the Pacific.
The direction of gene flow was estimated using a Bayesian coalescent approach implemented in MIGRATE 4.4.4 (Beerli et al., 2019), which infers the mutation-scaled effective population size θ (= N e µ for mtDNA data, where N e is effective population size and µ is mutation rate) and migration rate M (= m/µ,  Table 1. Different colors refer to different haplotypes. The landmass exposed during the Last Glacial Maximum (Ludt and Rocha, 2015) is shaded light gray. The directions of gene flow estimated between adjacent E. ferox populations are indicated with black arrows. Only six representative populations were examined (J: Japan; HPI: Hoping Island; SM: Shimen; DT: Datan; LM: Longmen; Dongyuping). Gene flows are presented by the effective number of migrants per generation (N e m), with 95% confidence intervals shown in parentheses.
Frontiers in Marine Science | www.frontiersin.org where m is the immigration rate per generation). Based on these parameters, the effective number of migrants per generation was then calculated as N e m = θM (Beerli et al., 2019). Because the distance between Malaysia/Singapore and the remaining populations is too far (>2,000 km Euclidean distance apart), samples of those two regions were excluded from gene flow estimation. Of the remaining populations, six were selected for analysis, which being two from Penghu (Dongyuping and Longmen), three from northwestern Taiwan (Datan, Shimen and Hoping Island), and one from Japan. To perform coalescent migration analyses, a HKY model was applied, and the other priors were set as default. The MCMC simulation consisted of one long chain with 500,000 recorded genealogies after a 10% burn-in, a sampling increment of 100, and a static heating scheme with four chains and temperatures {1.0, 1.5, 3.0, 10 6 }. To test the stepping stone hypothesis, the analyses were conducted with (1) a stepping stone model with symmetric migrations and (2) a migration matrix model specified to allow full migrations between adjacent populations.

Sequence Diversity and Population Differentiation
All 317 samples were amplified and sequenced successfully. The length of raw mtCOI sequences ranged from 499 bp to 536 bp. After trimming, 499 bp were retained for the following analyses. Nucleotide composition profiling demonstrated that sequences were A + T-rich (A, 33.43%; T, 27.38%). After alignment, 41 variable sites were detected, predominantly transition substitutions (Ti: Tv = 43.11), 16 of which were phylogenetically informative. A total number of 42 haplotypes were identified, with the highest found in Wangan ( Table 1). The relative abundance of the haplotypes in each site is shown in Figure 1. High haplotype diversity (h) was detected in all populations, with a mean value of 0.814 ± 0.014. In contrast, nucleotide diversity (π) was relatively low, with a mean value of 0.00337 ± 0.00014 (Table 1). Median-joining network (MJN) analysis revealed a star-like network, consisting of one major dominant haplotype, from which two second dominant haplotypes and numerous singletons derived. The latter two also contained several derived singletons of their own. All haplotypes differed from each other by only one or two mutation steps (Figure 2). No distinct haplotype clustering with respect to geographical location was detected. Pairwise F ST comparisons also suggested low genetic differentiation among E. ferox populations ( Table 2), indicating little or no phylogeographical structure in E. ferox.

Historical Demography and Gene Flow Estimation
Neutrality tests indicated that E. ferox experienced a sudden population expansion, with support from significant negative values of Tajima's D (−2.08071, p < 0.01) and Fu's Fs indices (−27.53360, p < 0.001), together with a significantly positive value of Ramo-Onsins and Roza's R 2 index (0.01973, p < 0.05) in a combined analysis of all populations in the northwestern Pacific (Table 3). At the regional scale, results showed that all FIGURE 2 | Median-joining network of Eriphia ferox. Each circle represents a unique haplotype. Circle sizes are proportional to haplotype frequencies shown in parentheses. Each slash across a line connecting haplotypes indicates a single mutation step.   Table 1.
Frontiers in Marine Science | www.frontiersin.org populations experienced sudden population growth (significant negative D and Fu's Fs, significant positive R 2 ) except for Malaysia (Table 3). Population expansion in E. ferox was further illustrated by unimodal mismatch distributions (Figure 3). Non-significant SSD statistics and rg indices were observed in populations from Taiwan, Malaysia, and Singapore, implying that the polymorphism pattern in these populations conformed to unimodal frequency distribution, except for Penghu (Table 3).
An SSD estimate and rg index were not available for the E. ferox population in Japan due to highly differentiated sequences that Arlequin failed to process. Bayesian skyline plot (BSP) analysis revealed that the initiation of overall E. ferox population growth occurred approximately 10,000 years ago, and the effective population size reached a plateau around 2,000 years ago (Figure 4A). At the regional scale, the effective population sizes of E. ferox remained stable in Japan, Malaysia, and Singapore, probably because the sample sizes of those populations were too small (≤20 individuals each) to estimate fluctuation through time (Grant, 2015). In contrast, populations in Taiwan and Penghu both experienced recent population expansions over the past 20,000 years. The Taiwanese population started expanding roughly 16,000 years ago, and reached a stable effective population size about 6,000 years ago ( Figure 4B). Population expansion in Penghu began approximately 10,000 years ago, and saturated ca. 2,000 years ago ( Figure 4C). For coalescent migration analyses, the migration matrix model allowing full migrations between adjacent populations (model 2) was selected as the more probable model over the stepping stone model (model 1) using log Bayes factors (LBFs). The log marginal likelihood for model 1 and model 2 were −1134.916 and −1113.828, respectively. The value of LBF was −42.176. Hence, the stepping stone hypothesis was rejected. Our results demonstrated a highly asymmetric gene flow between adjacent E. ferox populations along the east Taiwan Strait (Figure 1). Centered at Datan, gene flow occurred in two directions: one migrated northward with approximately 1.23 effective migrants (N e m) per generation from Datan to Shimen, 1.4 N e m from Shimen to Hoping Island, and then 1.5 N e m from Hoping Island to Japan (Figure 1). The other moved southward with 1.37 N e m from Datan to Longmen, and then 1.71 N e m from Longmen to Dongyuping in Penghu (Figure 1).

DISCUSSION
The Taiwan Strait plays an important role in shaping marine phylogeography across marginal seas in the northwestern Pacific. Yet, previous broad-scale research in the region preferentially emphasized the west Taiwan Strait (i.e., Chinese coast), leaving the east Taiwan Strait barely explored. In the present study, we analyzed 317 partial mtCOI sequences of the ferocious reef crab Eriphia ferox collected from 21 sites, spanning nearly 4,500 km (Euclidean distance) within the northwestern Pacific. With specific focus on the populations inhabiting the northwestern Taiwan and Penghu archipelago, we unveil the detailed historical connectivity and demography of E. ferox through the east Taiwan Strait. Despite mtCOI reveals only partial information, it is still the benchmark for non-model organisms lack of genetic information, e.g., E. ferox in this study. Our analyses of the haplotype network and pairwise F ST comparisons demonstrated a lack of phylogeographical structure in E. ferox populations, implying the existence of a migration corridor connecting the South and East China Seas (SCS and ECS) through the east Taiwan Strait. Previous studies along the Chinese coast of the west Taiwan Strait have reported significant genetic structures, e.g., mollusks and crustaceans (reviewed in Ni et al., 2014), suggesting that there is no such corridor along the west Taiwan Strait. Our discovery of the corridor along the east Taiwan Strait is probably the first ever reported in this region. Moreover, instead of serving as a stepping stone, our analysis of coalescent migration demonstrated that Taoyuan Algal Reef, the only biotic reef in the west coast of Taiwan, is not just essential, but also the only population source that exports migrants both northward and southward to adjacent populations. Our results indicated that northward gene flow propagated sequentially from Datan to Shimesn, Hoping Island, and then Japan. Southward gene flow extended from Datan to Longmen (northern Penghu) and then Donyuping (southern Penghu). Previous studies indicated that the migration rates approximating one effective migrant per generation that we estimated in this study were sufficient to maintain genetic diversity and population viability (Wang, 2004;Nathan et al., 2017). The migration of benthic invertebrates mostly depends on larval dispersal via ocean currents. The panmictic population structure of E. ferox across the northwestern Pacific observed in our study suggested its potential for long distance dispersal. Although so far there is no estimate of the pelagic larval duration for E. ferox, larvae of its congeneric species E. verrucosa (Forskål, 1775) were reported to last for 49 days, contributing to remarkably long-distance connectivity across the East Atlantic and Mediterranean Sea (Deli et al., 2019).
Based on the information above, we speculated that the lack of migration corridor along the west Taiwan Strait results from the difference in coastal hydrology along the two sides of the Taiwan Strait while the dispersal pattern of E. ferox centered at Datan is shaped by the complex current system in the east Taiwan Strait. In the east Taiwan Strait, the Kuroshio's branch or loop current flows unidirectionally northward year-round along the west coast of Taiwan through the Penghu channel (Hu et al., 2010). Both seawater temperature and salinity are  rather warm and stable in the east Taiwan Strait, providing a smooth pathway for the dispersal of tropical/subtropical marine faunas like E. ferox. The reproduction season of E. ferox is reported to occur from July to October (previously recognized as E. smithii in Tomikawa and Watanabe, 1992), thus the northward gene flow from Datan should be attributed to the conveyance of pelagic larvae by Kuroshio's loop current along the east Taiwan Strait in summer. Conversely, the southward gene flow originating from Datan is thus more likely to arise from secondary contact due to post-LGM population expansion. The dispersal direction from Taiwan to Penghu is congruent with our BSP estimation that population expansion in Taiwan occurred before the population expansion in Penghu. Contrarily, surface currents in the west Taiwan Strait are influenced more substantially by regional monsoons, which reverse in direction on a semiannual basis. In summer, the extension of the South China Sea Warm Current (SCSWC), together with the Yuedong Coastal Current driven by the southwestern monsoon, flows northward along the Chinese coast in the west Taiwan Strait (Hu et al., 2010). In winter, the counter-wind SCSWC tends to be weaker or even disappears in the upper layer and current direction changes from northward to southwestward when the southwestward Zhemin Coastal Current (ZCC) propelled by the northeastern monsoon dominates in the region. Moreover, the ZCC brings in diluted cold water originating from the Changjiang estuary (Hu et al., 2010), creating a steep sea surface temperature (SST) gradient, which in turn forms a dispersal barrier in the west Taiwan Strait and contributes to the genetic differentiation of lineages between the SCS and ECS along the Chinese coast (Ni et al., 2012;Xu and Chu, 2012;Wang et al., 2017;Qu et al., 2021).
We elucidated a rapid population growth following a population bottleneck throughout E. ferox demographic history, inferring from the combination of high mean haplotype diversity and low mean nucleotide diversity (Grant and Bowen, 1998;Avise, 2000). The recent or sudden population growth of E. ferox is supported by our demographic analyses of neutrality test and mismatch distribution. Bayesian skyline plots (BSPs) further indicated that the E. ferox population expanded during the post-LGM period, when a rapid rise in sea level submerged the Taiwan Strait and stabilized to the present level approximately 7,000 years ago (Kimura, 2000;Li et al., 2014). Such population expansion of E. ferox occurred first at Taiwan and followed at Penghu, possibly reflecting the direction of seawater intrusion into the east Taiwan Strait after the end of the LGM. The demographic expansion of E. ferox also coincided with the time that Taoyuan Algal Reef started to form around 7,500 years ago. Hence, E. ferox might have colonized Taoyuan Algal Reef since then and gradually grew to its present substantial population size.
In conclusion, our phylogeographical study of E. ferox based on mtCOI sequences identifies a potential migration corridor connecting SCS and ECS through the east Taiwan Strait. Instead of serving as a stepping stone, we showed that Taoyuan Algal Reef, represented by Datan Algal Reef, is confirmed as an essential and also the only population source along this migration pathway for reef-associated biotic species in the northwestern Pacific. To complete the picture of marine phylogeography in this region, it is necessary to involve more reef-associated marine taxa with different life histories in future studies. The gaps between Penghu and Malaysian/Singapore populations should be filled by increasing the number of sampling sites in between. The pattern of demographic history and genetic connectivity of E. ferox described in this study should also be compared by additional nuclear markers and multi-locus genetic approaches; e.g., microsatellites or single nucleotide polymorphisms (SNP). Still, our study, along with other pioneer research (Kuo et al., 2020;Heard et al., 2021), demonstrates the uniqueness and pivotal role of the Taoyuan Algal Reef along the east Taiwan Strait. Such a distinctive ecosystem with an ancient origin is now confronting catastrophic disturbance. The construction of an industrial port and terminals to receive and store liquefied natural gas has undermined the integrity of the algal reefs, leaving the survival of its marine inhabitants in jeopardy (Liou et al., 2017;Silver, 2018). As this ecosystem is severely understudied to date, further research is urgently needed to reduce the impacts of construction and enable conservation management strategies that ensure the sustainability of Taoyuan Algal Reef.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/genbank/, OK509858-OK510174.