Causations of phylogeographic barrier of some rocky shore species along the Chinese coastline

Substrate, ocean current and freshwater discharge are recognized as important factors that control the larval dispersal and recruitment of intertidal species. Life history traits of individual species will determine the differential responses to these physical factors, and hence resulting in contrasting phylogeography across the same biogeographic barrier. To determine how these factors affect genetic structure of rocky shore species along the China coast, a comparative phylogeographic study of four intertidal and subtidal species was conducted using mitochondrial and nuclear DNA by combining new sequences from Siphonaria japonica with previously published sequences from three species (Cellana toreuma, Sargassum horneri and Atrina pectinata). Analysis of molecular variance and pairwise ΦST revealed significant genetic differences between the Yellow Sea (YS) and the other two marginal seas (East China Sea, ECS and South China Sea, SCS) for rocky-shore species (S. japonica, C. toreuma, S. horneri), but not for muddy-shore species Atrina pectinata. Demographic history analysis proved that the population size of all these four species were persistent though the Last Glacial Maximum (LGM, ~20 ka BP). Migration analysis revealed that gene flow differentiated northward and southward migration for these four species. However, the inferred direction of gene flow using alternatively mitochondrial or nuclear markers was contradictory in S. japonica. It is concluded that there is a phylogeographical break at the Yangtze River estuary for the rocky shore species and the causation of the barrier is mainly due to the unsuitable substratum and freshwater discharge. All four intertidal and subtidal species appear to have persisted through the LGM in China, indicating the lower impact of LGM on intertidal and subtidal species than generally anticipated. The imbalanced gene flow between YS and ESCS groups for these four species could be explained by historical refugia. The discordance between mitochondrial and nuclear markers in the MIGRATE analysis of S. japonica prove the importance of employing multi-locus data in biogeographic study. Climate change, land reclamation and dam construction, which are changing substrate and hydrological conditions around Yangtze River estuary, will consequently affect the biogeographic pattern of intertidal species.


Background
Phylogeographical patterns of marine faunas are complex and affected by multiple biotic and abiotic factors, and it is in the long-term interest of marine ecologists to understand the roles of these factors in determining the distribution and genetic structuring of species. Glacialinterglacial climate fluctuations during the Pleistocene led to changes of sea level [1,2] and then caused habitat contractions or expansions [3], which impeded gene flow of marine species and resulted in genetic divergence. Following the glacial retreat, demographic expansion occurred in most marine taxa [4]. Postglacial exchanges of propagules may erase the signals of historic isolation. Accordingly, the biology of the species (e.g. dispersal capacity), availability of suitable habitat and ocean current regimes, which determine the contemporary level of gene flow, would significantly contribute to promote, maintain or homogenize the genetic divergence created by glacial periods. Suitable habitat is assumed to have a significant impact on the phylogeographical distribution of coastal species [5][6][7][8][9][10] .
The impacts of ocean currents on larval dispersal are variable. Sometimes oceanic currents can promote larval dispersal and enrich population connectivity [11][12][13][14]. However, converging ocean currents may also pose a potential barrier to gene flow to some extents [15]. Additionally, a river outflow carrying large amount of sediments and freshwater discharge can influence physical and chemical characteristics (e.g. geomorphology, turbidity, salinity, nutrients and dissolved oxygen etc.) of estuarine and coastal water, which control the biotic dynamics in estuary and coastal ecosystems [16]. This has significant impacts on the dispersal of marine taxa [5,15,[17][18][19][20].
The marginal seas in the northwestern Pacific have changed dramatically in area and configuration, particularly during the Pleistocene glacial-interglacial cycles [21]. Three of the marginal seas, namely the South China Sea (SCS), the Yellow Sea (YS) plus the East China Sea (ECS), and the Sea of Japan (SOJ), served as three independent glacial refugia and resulted in vicariant divergence in marine fauna [22][23][24]. On the other hand, the lack of genetic structuring in some sympatrically distributed taxa was attributed to their postglacial colonization of these regions [24]. Yet, these hypotheses were not rigorously tested using comprehensive phylogenetic data because most of the studies performed in the region were based on a single marker (in most case the mitochondrial DNA) and largely biased to commercial species, of which the effect of anthropogenic introduction for aquaculture purposes was unknown. Hence, additional investigation using both mitochondrial and nuclear markers was advocated [24].
The Yangtze River represents another critical factor in the gene flow of coastal species in China. The Yangtze River, the fifth largest river in the world in terms of volume discharge, brings huge amounts of water (8~9 × 10 11 m 3 ) into the East China Sea annually [25]. When the freshwater discharge reaches its maximum in spring and summer, Yangtze River plume can extend eastward or northeastward to the Jeju Island (126°08′~126°5 8′E, 33°08′~33°60′N) [26,27], and it can dramatically change surrounding ocean currents [26,28] and salinity of the upper layer of the Kuroshio Current [29,30] (Fig. 1). Moreover, sediment discharges from the Yangtze and other nearby rivers have formed the Yangtze River Delta with an area of more than 3 × 10 4 km 2 [31,32]. The~600 km long coastline from Lianyungang, Jiangsu Province (34°36′N, 119°13′E) to Shaoxing, Zhejiang Province (30°19′N, 120°4 6′E) is mainly salt marsh [33]. This together with the Yangtze River discharge is assumed to form a contemporary dispersal barrier for marine species that require a hard substratum (e.g. intertidal rocky shore) or with larvae that cannot tolerant decreased salinity, including the gastropod Cellana toreuma [20], the bivalve Cyclina sinensis [34] and the macroalga Sargassum hemiphyllum [18,35]. However, studies on other coastal fauna did not detect any genetic split across the Yangtze River [24]. Therefore, more comprehensive study is required to test for the influence of this barrier.
To detect causations (such as substrate, freshwater discharge, coastal current and historical events) of phylogeographic break of intertidal and subtidal rocky shore communities, a multi-species genetic analysis of several intertidal and subtidal species along the Chinese coast was performed. We combined new sequences with  [30] previously published sequences from the Chinese coast to construct a database of mitochondrial and nuclear DNA from four species, including three rocky-shore species (limpet Siphonaria japonica, limpet C. toreuma [20], macroalga S. horneri [36]) and one muddy-shore species (bivalve Atrina pectinata [37]). These species possess different habitat preferences, reproductive seasons and breeding modes (See Table 1, [37][38][39][40]), which provide insights into the responses of intertidal and subtidal species to multi-factors.

Sequence variations Siphonaria japonica
High levels of haplotype diversity (h, mean ± S.D.) could be observed in all nine locations of Siphonaria japonica for both markers (COI, 0.902 ± 0.052~0.989 ± 0.019; ITS, 0.986 ± 0.016~1.000 ± 0.009; Table 2). The difference in haplotype diversity between four locations from the Yellow Sea (YS) and five locations from the East China Sea (ECS) and the South China Sea (SCS) was statistically significant in COI gene (P < 0.05), but not in ITS sequences.

Cellana toreuma, Sargassum horneri and Atrina pectinata
Dong et al. [20] had suggested that the two populations of C. toreuma from YS have significant higher haplotype and nucleotide diversity as compared to the populations from ECS and SCS. Re-analysis of published data of S. horneri revealed that the difference between populations from YS and populations from ECS/SCS was significant in nucleotide diversity (P = 0.034), but not in haplotype diversity (P = 0.187). For A. pectinata, however, there were no significant difference between populations from YS and populations from ECS in both haplotype diversity and nucleotide diversity.

Population structure Siphonaria japonica
The global Φ ST analysis of mtDNA and nuclear DNA in S. japonica exhibited values of Φ ST significantly different from zero (Additional file 1: Table S1), indicating a significant spatial genetic structure. Pairwise Φ ST analysis of COI and ITS sequences among locations showed that significant genetic differentiation (P < 0.05) between locations from the YS and locations from the ECS and SCS (Table 3). Pairwise Φ ST values were also calculated among different groups. The YS group was significantly different from the ECS group and the SCS group, and the difference between the ECS group and the SCS group was low and non-significant (Table 4), which is similar to the results obtained from C. toreuma [20].
We tested the hypothesis of reduced gene flow between YS and ECS/SCS as observed in pairwise Φ ST . Therefore, the locations of the four species were subdivided into two groups (YS group including locations from Yellow Sea and ESCS group containing locations from East and South China Seas) for AMOVA analyses. The hierarchical analysis of AMOVA for S. japonica indicated that variation within locations (Φ ST ) accounted for 70.09 % (P < 0.001) and 93.10 % (P < 0.001) for COI and ITS respectively, followed by variation among groups (Φ CT ) (CO1: 30.27 %, P = 0.005; ITS: 6.87 %, P = 0.007) ( Table 5). The variations among locations within groups were −0.36 % (P = 0.827) for COI and 0.03 % (P = 0.422) for ITS (Table 5).

Cellana toreuma, Sargassum horneri and Atrina pectinata
One-group AMOVA indicated that the genetic variation among all samples was negative and insignificant in A. pectinata (Φ ST = −0.01), but positive and significant in C. toreuma and S. horneri (Additional file 1: Table S1). At location level, re-analysis of pairwise Φ ST showed that there were no significant differences between all eight locations in A. pectinata (Additional file 2: Table S2). However, there existed significant difference between locations from YS and locations from the ECS and SCS (except for FJ location) (Additional file 3: Table S3) in S. horneri. At group level, the YS group was significantly different from the ECS group and the SCS group and the difference between the ECS group and the SCS group was low and non-significant in S. horneri (Table 4). However, there was no significant difference between the YS group and the ECS group for A. pectinata (Φ ST = −0.0001, P = 0.436, Table 4).
Under the grouping criteria mentioned above, AMOVA analysis for COI of C. toreuma showed that there were significant genetic variations among groups (19.39 %, P < 0.001) and within locations (80.82 %, P < 0.001) ( Table 5). AMOVA analysis for COIII of S. horneri showed significantly high variations among locations within groups (15.14 %, P < 0.001) and variation within locations (85.49 %, P < 0.001). The variance among groups was relatively low and insignificant (0.37 %, P = 0.331; Table 5). For COI sequences of A. pectinata, variance components at all three levels were non-significant (Table 5).

Phylogenetic analysis Siphonaria japonica
Though branches related to geography were not significantly deep, two putative groups were suggested in the median-joining network and in neighbour-joining tree of mtDNA COI of S. japonica ( Fig. 2a; Additional file 4: Figure S1). One group (northern group) contained all individuals except one (WH22) from the Yellow Sea and 33 individuals from the East and South China Seas. The second group (southern group) included the majority of East and South China Seas individuals plus the single specimen from WH (WH22). Based on network and NJ tree, haplotype relationships of ITS sequences revealed no significant branches or clusters corresponding to geography ( Fig. 2b; Additional file 5: Figure S2).

Cellana toreuma, Sargassum horneri and Atrina pectinata
The haplotype networks of other three species (C. toreuma, S. horneri and A. pectinata) showed a pattern which did not exhibit obvious subdivision according to geographical locations (Fig. 2c, d and e).

Demographic analysis Siphonaria japonica
Mismatch distribution analysis showed a unimodal distribution for both northern group and southern group of COI in S. japonica, suggesting rapid demographic   Fig. 3a and b). The tau value (τ) can provide a rough estimation of the time when rapid population expansion began. Using a substitution rate of 1 % per million years [9,41], the time of demographic expansion for northern group and southern group was estimated as~440 ka (90 % CI: 313.9-682.4 ka) BP and~190 ka (90 % CI: 104.4-259.3 ka) BP, respectively ( Table 6). The BSPs indicated that the northern group had a slightly increase of population size since 475 ka BP, but the most distinct demographic expansion happened around 200 ka BP (Fig. 4a). The southern group exhibited a gentle growth in population size from~225 ka BP (Fig. 4b), when the earth was at an interglacial stage (Fig. 4c, [42,43]). The inconsistent estimates of demographic expansion could be the difference between modeling an expansion with pairwise differences in the mismatch distribution analysis versus the coalescent theory in the BSP analysis. A posterior simulation-based analogue of Akaike's information reiteration through MCMC (AICM) test for northern group revealed that the expansion model favored over the other two models (Additional file 6: Table S4), and the BSP showed a demographic expansion. The constant size model was the best-fitting demographic model for southern group (Additional file 6: Table S4) which exhibited a similar pattern to the BSP.

Cellana toreuma, Sargassum horneri and Atrina pectinata
Although mitochondrial sequences in both C. toreuma and S. horneri showed a significant value of Φ ST , we reconstructed their demographic history with mismatch distribution and BSP because, no spatial subdivision of haplotypes of mitochondrial sequence existed in the network (Fig. 2c, d). The BSP analysis showed that the population size of C. toreuma had an extremely gentle increase since 260 ka BP (Additional file 7: Figure S3A), and the population size of S. horneri experienced a very slow increase since 194 ka BP with the mutation rate of 1.675 % MY −1 (Additional file 7: Figure S3B). The population size of A. pectinata increased significantly since 232 ka BP with the mutation rate of 0.775 % MY −1 (Additional file 7: Figure S3C). Mismatch distribution analysis showed that the shapes of mismatch distribution  [20]; Sargassum horneri, Hu et al. [36]; Atrina pectinata, Liu et al. [37].*P < 0.05; ***P < 0.001 Rocky    were unimodal (Fig. 3c, (Table 6). Therefore, the demographic expansion of the two groups of S. japonica and other three species significantly predated the Last Glacial Maxima (~20 ka BP).
AICM tests for C. toreuma and S. horneri showed that the constant size model was a better fit to the data than the BSP and expansion models (Additional file 6: Table S4). It was not surprising for these two rockyshore specie as they both showed relatively flat BSPs.
The expansion model was strongly favored for A. pectinata that showed a significant demographic expansion in BSP.

Gene flow Siphonaria japonica
Because the present study mainly focused on the gene flow across the Yangtze River Estuary, locations of each species were divided into YS group and ESCS group to detect the possible phylogeographical break. In the present studies of S. japonica, contrasting patters of gene flow between mitochondrial and nuclear markers were suggested ( Fig. 2a and b, Additional file 8: Table S5). For

Discussion
Phylogeography in intertidal and subtidal species along the China coast: northern and southern differentiation Significant population differentiation existed among northern (Yellow Sea) and southern groups (East China Sea and South China Sea) of Siphonaria japonica, Cellana toreuma and Sargassum horneri respectively, indicating that there is a genetic break between northern and southern populations of rocky-shore species. For mitochondrial and nuclear sequences of S. japonica, pairwise Φ ST analyses among different locations indicated that Yellow Sea locations were significantly different from locations from the East China Sea and the South China Sea (Table 3). When these locations were divided into three groups (YS, ECS and SCS) based on their geographic locations, the Φ ST values between the YS group and the other two groups (ECS and SCS) were significant for limpet S. japonica and macroalga S. horneri. This result is similar to results for the limpet C. toreuma, and suggests that the Yellow Sea group is relatively isolated from the other two groups. AMOVA analyses for COI and ITS of S. japonica and for COI of C. toreuma showed that genetic differentiation among YS and ESCS groups accounted for a high proportion of the molecular variance among groups (Table 5), which suggested that populations were mainly grouped according to the Yangtze River Estuary. Phylogenetic analysis showed that two putative groups existed in COI of S. japonica, which also demonstrated a genetic break around the Yangtze River Estuary. A geographic genetic break can also be detected by changes in gene diversity between populations. Significant difference in genetic diversity between locations from YS and locations from ECS/SCS for these three rocky-shore species supported the suggestion that a genetic break existed between YS and ECS/SCS groups. In contrast to the rocky-shore species, no obvious population structure for the muddy-shore species (Atrina pectinata) was found based on the results from AMOVA analysis and pairwise Φ ST values. The absence of genetic difference between YS and ECS groups as detected for the muddy-shore species was similar to previous studies of Rapana venosa [44], Tegillarca granosa [45] and Cyclina sinensis [46].
During the Pleistocene glacial-interglacial cycles, areas and configurations of the marginal seas have changed dramatically in the West Pacific [21,47]. When the sea level fell about 120-140 m over the past~800 kyr [48], the Yellow Sea (YS) and the East China Sea (ECS) were reduced to an elongated trough, the Okinawa Trough, and the South China Sea (SCS) became a semi-enclosed Table 6 Results of population expansion tests on the four species. Fu's F S , Tajima's D and mismatch distribution estimate (τ) and the real expansion time (t) with 90 % credibility intervals in parentheses for each species are included. All COI sequences of S. japonica were divided into two groups: northern group and southern group, according to the phylogenetic analyses (See the results in the phylogenetic analysis). Reference: Cellana toreuma, Dong et al. [20]; Sargassum horneri, Hu et al. [36]; Atrina pectinata, Liu et al. [37]. *P < 0.05; **P < 0.01, ***P < 0.001  gulf [49]. ECS and SCS were separated by a land bridge which connected Taiwan and the continent [50]. Previous comparative phylogeographical studies in marginal seas of the northwestern Pacific have suggested that historical isolation between SCS and ECS plays important roles for the present-day distribution of genetic variation of coastal species such as some fishes, crustaceans and muddy/sandy shore mollusks [24].
Multi-factors controlling the population genetic differentiation of rocky-shore species The contrasting phylogeographical patterns between rocky and muddy intertidal species indicates that substrate plays important roles in the phylogeographical patterns of intertidal species. From Lianyungang, Jiangsu Province to Shaoxing, Zhejiang Province, the absence of appropriate habitats (~600 km salt marsh shore) could hamper the settlement of rocky intertidal species and consequently genetic exchange between Yellow Sea and East China Sea populations. The negative impacts of a lack of appropriate substrate on genetic connectivity have been reported in other intertidal fauna. For instance, the existence of long stretches of sandy beach serves as a barrier to dispersal between Cape St Lucia and Zinkwazi Beach in the limpet S. nigerrima in southeast Africa [9]. Recently, a genetic analysis of S. japonica specimens collected from Yangguang Island (32°36′N, 121°08′E), an artificial island between Zhoushan (30°01′N, 122°06′E) and Lianyungang (34°36′N, 119°13′E), was carried out and revealed that this location is genetically similar to the location in Zhoushan. These results indicate that colonization of some rocky intertidal species can happen across the barrier if suitable habitat is provided (Huang XW, Wang W and Dong YW, unpublished observations). In contrast, as seen above, the species living in muddy substrate appear unaffected by the salt marsh around at the Yangtze River estuary. Freshwater discharge influences hydrological condition nearby estuary and may have some impacts on gene flow between the Yellow Sea and East China Sea. During spring and early summer, the spawning season of S. japonica, the Taiwan Warm Current (TWC) and the China Coast Current (CCC) in ECS and SCS flowing northward [51] transport pelagic larvae from SCS to ECS. However, the freshwater discharge from the Yangtze River will cause deflection of the East China Sea coastal current. The size and distance of deflection is more prominent in spring and summer due to the increased amount of surface runoff in the rainy season [52]. In spring and summer, the plume of water from the Yangtze River discharge shifts in a northerly direction in parallel with the TWC with a clockwise deflection (Fig. 1b). When it reaches its maximum flow, the Yangtze River discharge affects surrounding hydrological conditions [26,28] and causes a decline in the salinity of the upper layer of the Kuroshio Current [29], which could influence the northward transport of larvae into the Yellow Sea. Dong et al. [20] also suggested that unique haplotype and higher genetic diversity in YS group were mainly contributed to the ocean current and freshwater discharge during the spawning season of C. toreuma. The impact of freshwater discharge on the phylogenetic distribution of marine species has been widely observed in previous studies. For example, the outflow from the Amazon River has been invoked as a major factor for the biogeographic break between Brazilian and Caribbean faunas (e.g. barnacle Chthamalus proteus [53] and surgeonfishes Acanthuridae sp [5]).
Although the same phylogeographic break around the Yangtze River estuary was observed for two rocky intertidal limpets (C. toreuma [20]; S. japonica, the present study), haplotype networks of these two limpets were different. The haplotype network of C. toreuma is a single star-like network [20] and S. japonica presents a relatively complex network with two putative groups (Fig. 2). C. toreuma lays its eggs into sea water directly and the length of the pelagic larval stage varies between 4-18 days in congenerics [54]. S. japonica deposits gelatinous egg ribbons containing numerous small eggs on the rocky shore, and then the eggs develop into planktotrophic veliger larvae within~15 days. The veliger larvae of S. japonica could at least maintain in the plankton about seven days (Wang W and Dong YW, unpublished observations). Therefore, the discrepancy of phylogeographic patterns between S. japonica and C. toreuma could be associated with their different reproductive modes and larval dispersal capabilities. The high dispersal capability of the brown macroalga S. horneri could be partly attributed to low variance among groups based on the AMOVA analysis. S. horneri can breed by both sexual reproduction and asexual reproduction, and can form floating mats which can drift about 1-5 months after being detached from the substratum [55][56][57].

Historical demography of intertidal and subtidal species along the China coast
Even though the expansion model and the constant size model received strong support over the BSP for the putative northern and southern groups in S. japonica, respectively, evidences from the mismatch distribution analysis and the BSP suggested northern group and southern group had different timings of population expansion events. The expansion time of northern group (313.9~682.4 ka BP) was earlier than that of southern group (104.4~259.3 ka BP). This corroborates recent evidence of the existence of a northern refugium in the Northwestern Pacific observed in other marine organisms (e.g. seaweed Ishige okamurae [58], limpet C. toreuma [20] and barnacle Chthamalus challengeri [59]). In northwestern Pacific, South China Sea (SCS) and East China Sea (ECS) are widely accepted as southern glacial refugia [22][23][24]. Thus, new evidences of the existence of a northern refugium in this area may be helpful to illustrate the influence both of past populations in glacial refugia and of contemporary gene flow in shaping current phylogeographical patterns in the future studies. Furthermore, the population sizes of these four species were persistent through the Last Glacial Maximum (LGM;~20 ka BP). Recent metaanalyses of demographic history of intertidal rocky organisms in the northeastern Pacific [3] and northwestern Pacific [24] also converged on similar finding that majority of the species were not extirpated entirely by the LGM and regional persistence maybe more prevalent.
At the contemporary level, MIGRATE analyses revealed that migration rates were unbalanced between ESCS and YS groups of S. japonica. However, the inferred direction of gene flow using alternatively mitochondrial or nuclear markers was contradictory. COI analyses suggested that southward migration rate was almost ten times the northward migration rate. ITS sequences suggested the northward migration rate was significantly higher (Fig. 2b). Discordances between nuclear and mitochondrial data in animal biogeographic studies are not uncommonly reported [60]. Adaptive introgression, demographic disparity and sex biased migration are commonly invoked as potential explanation of discordance observed [60]. Sex biased gene flow appears unlikely for this hermaphroditic limpet with planktonic larvae transported by ocean current while adaptive introgression is difficult to test based on the present data. The mitochondrial genes have a higher mutation rates and smaller effective population size than the nuclear genes, and hence considered to be more informative in shallow relationships. Besides, the high imbalance between immigration and emigration rates in other three species (Fig. 2c, d and e) could be associated with historical refugia and suggested the likely existence of sources and sinks at metapopulation level. Recent studies emphasize the need of multi-locus data for accurate estimation of various population parameters [3,61]. As a result, more data are required to determine the underlying mechanism leading to the discordance, and our results emphasize the importance of employing multi-locus data in biogeographic study.

Biogeography of intertidal species in China
The phylogeographic barrier of rocky intertidal species around Yangtze River estuary will possibly disappear under the coupled impacts of climate change and human activities. Firstly, northward shift of rocky intertidal species is ubiquitous in the scenario of climate change [62][63][64][65]. The coastal sea surface temperature (SST) and extreme hot days have continued to increase from 1982 to 2010 [66,67]. The increasing temperature will potentially force the northward shift of intertidal species along China coast; Secondly, land reclamation has resulted in about 55 % loss of coastal wetland in China from 1949 to 2002 [68], and large numbers of artificial structures provide hard substrates in areas where these are generally absent and act as stepping stones to connect the Yellow Sea populations and East China Sea/South China Sea populations; Finally, about 50,000 dams constructed in the Yangtze River catchment from the late 1950 to 2003 have a storage capacity of 22 % of the annual water discharge (200 × 10 9 m 3 in 2003) and result in a strong decrease of sediment and freshwater discharge in spring and summer [69]. The changing hydrological condition around the Yangtze River estuary could enhance the possibility of larval dispersal across it. Overall, the potential northward shift of rocky intertidal species will change the biogeographic pattern along the Chinese coast and adaptive management should be considered for future management of rocky intertidal ecosystem in China.

Conclusions
A significant phylogeographic break, occurring around the Yangtze River Estuary, was observed for populations of rocky-shore species along the China coast. Substrate, ocean current and freshwater discharge are suggested as major factors determining the contemporary structure of rocky-shore species by limiting north-south dispersal of planktonic larvae. In addition, historical events and life history characteristics can also influence the contemporary phylogeographic distribution of intertidal and subtidal species. However, human activities are changing the habitat of intertidal species. Large-scale land reclamation activities can change the sedimentary substrate, and so provide suitable hard substrates for colonization by rocky-shore species. On the other hand, numerous dams constructed in the Yangtze River catchment decrease the riverine sediment supply to the sea and will impact on the environment of the Yangtze Delta and the nearby coastal ocean.

Sampling and sequencing
A DNA sequence database consisted of new sequences (cytochrome c oxidase subunit I gene, COI; internal transcribed spacer, ITS) from Siphonaria japonica (specimens were collected from nine rocky shore localities along the Chinese coastline between May 2012 and January 2013) and previously published data (mitochondrial gene) from three other intertidal and subtidal species in the Pacific Northwest (Fig. 2).
Partial sequences of COI and ITS in S. japonica were amplified by polymerase chain reaction (PCR). COI sequences were amplified with universal primers LCO1490 and HCO2198 [70] and ITS sequences were amplified using primers its-1d and its-4r [71]. PCRs were conducted in a 25-μL reaction volume containing 2.5 μL of 10 × buffer (Mg 2+ Plus), 2 μL of 2.5 mM dNTPs, 1 μL of each 10 mM primers, 0.25 μL (1.25 U) of Taq DNA polymerase and 200 ng DNA template. Amplification was initiated with denaturing at 95°C for 3 min, followed by 35 cycles of 95°C for 1 min, annealing at 40°C for COI and 54°C for ITS for 1 min and 72°C for 1 min and then a final extension at 72°C for 10 min. After visualizing the target amplicon in 1.5 % agarose gels, PCR products of COI gene were sent to a commercial company for sequencing (Invitrogen Biotechnology Co., Ltd., Shanghai, China). The ITS products were purified using the PCR Purification Kit (Aidlab Biotechnologies Co., Ltd., Beijing, China), ligated to pMD19-T Vector (TaKaRa Biotechnology, Dalian, China) and then transformed into competent cell of Escherichia coli DH5α (TaKaRa Biotechnology, Dalian, China). Finally, a positive clone per individual was sequenced in both directions using M13 and M17 primers by Invitrogen Biotechnology Co., Ltd. (Shanghai China).

Sequence variation and population genetic analysis
All sequences in Siphonaria japonica were edited by comparing both strands using DNAMAN 7 software (LynnonBioSoft, Quebec, Canada), and then aligned with MUSCLE [72] using MEGA 5 [73] with default settings. Standard molecular diversity indices including haplotype diversity (h) and nucleotide diversity (π), and neutrality test including Fu's F s [74] and Tajima's D [75] were calculated using ARLEQUIN 3.5 [76]. jModelTest 2.1.1 [77] was utilized to estimate the bestfitting substitution model and substitution parameters with the Bayesian information criterion (BIC) for each species (Table 3). For sequences from S. japonica, phylogenetic analyses were performed based on the neighbour-joining (NJ) approach using MEGA incorporating the Tamura-Nei model (TrN) [78] with corresponding gamma correction for COI and ITS (the closest model in MEGA to TVM + G is the TrN + G model for NJ tree building). 1 000 bootstrap replicates were carried out to assess the clade credibility of the NJ phylogram. A haplotype network illustrating genealogical relationships between haplotypes was constructed for each species using the median-joining (MJ) method with Network 4.6 [79].
Pairwise Φ ST measures were calculated to evaluate the levels of genetic differentiation and significance was estimated with 10 000 permutations using ARLEQUIN. Furthermore, hierarchical analysis of molecular variance (AMOVA) [80] was performed to test for possible phylogeographic separation. In the present study, genetic differentiation among locations was only analyzed in S. japonica, because the published data from other three species have been analyzed. Because this study focused on the gene flow and phylogeographic break, the locations of each species were divided into three groups: the Yellow Sea (YS) group, the East China Sea (ECS) group and the South China Sea (SCS) group based on the geographical locations (refer to Table 4 for the allocation of sites). Because no significantly genetic differentiation was observed between ECS and SCS (see the Results Table 3), the ECS and SCS were combined as the ESCS group to compare with the YS group to test for the hypothesis of reduced gene flow across the Yangtze River outflow. On the other hand, samples from all locations for each species were considered as a single group to verify the significance of partitioning of genetic variance among all samples.

Historical demography
The coalescent-based approach is widely used to reconstruct the demographic history. In the present study, the Bayesian skyline plot (BSP) was generated in BEAST v1.7.4 [81] to estimate the change in effective population size over time. The analyses were only performed with the mitochondrial dataset because no distinct clade was revealed in ITS and the mutation rate is unknown for the marker, hampering the estimate of temporal scale. Analyses were conducted under corresponding model suggested by jModelTest for each species with constant Bayesian skyline tree priors with 10 groups under a strict clock model. Default priors were used for parameter settings. Three independent MCMCMC searches were run for 100 million generations and parameters were recorded every 10 000 generations with the first 10 million generations discarded as burn-in. For each run, the effective sample sizes (ESS) of important parameters sampled from the MCMCMC were >1000 in all three replicate runs calculated by TRACER 1.4 [81]. Results from three runs were then pooled with LOGCOMBINER version 1.7 [81] for final reconstruction of the BSP. Pairwise mismatch distribution analysis was also performed using the mitochondrial dataset of each species with ARLEQUIN 3.5 [76] to validate the result from BSP. Population expansion parameters (τ) were transformed to estimates of realtime since expansion (t) with the formula τ = 2μkt, where μ is the mutation rate and k is the sequence length. As two geographically-partitioned groups were revealed in phylogenetic analysis of COI in S. japonica (Fig. 2a), the demographic history of each putative group was reconstructed using both methods above.
Because of the absence of clear fossil or geological record, Colgan & Costa [41] and Teske et al. [9] had used a mutation rate of 1 % per million years to estimate divergence times in Siphonaria genus, when referring to available calibrated fossil data for marine gastropod.
Thus, in this study, such a mutation rate (1 % Myr −1 ) with a generation of 1 year was applied for COI of S. japonica. A divergence rate of 0.85-1.15 % Myr −1 for COI was set in C. toreuma, as used for estimation times of the C. nigrolineata [8]. The divergence rate was set at 2.6-4.1 % Myr −1 for COIII of S. horneri [36] and 0.7-2.4 % Myr −1 for COI of A. pectinata [37], respectively, which were proposed in the original papers. In addition a generation of 1 year for C. toreuma [39] and S. horneri [82] and a generation of 2 years for A. pectinata [83] were assumed. In the present study, the mutation rate for these three species (C. toreuma, 0.5 % Myr −1 ; S. horneri, 1.675 % Myr −1 ; A. pectinata, 0.775 % Myr −1 ) was obtained by averaging the divergence rates used above and then dividing by two.
To evaluate whether the BSP was the best model in reconstructing demographic histories, BSP model was compared with two simple moles in Tracer V 1.6: constant population size and expansion growth. The two alternative models were run in BEAST using the same way as mentioned above for BSP. A posterior simulation-based analogue of Akaike's information criterion through MCMC (AICM) [84] was used to compare all three models, which measured AIC from the posterior of each model, with score > 10 as strong evidence in favor of one model over the others [85].

Gene flow analysis
To examine the level and direction of contemporary gene flow across the phylogeographic break, we adopted the coalescent-based approach implemented in MIGRATE-N 3.5.1 [61] for estimating the migration rates between groups within the four studied species. Locations were grouped into a "YS" group and an "ESCS" group according to their geographical distribution. Random sub-samples were performed to allow the two groups to contain comparable number of individuals. The Bayesian approach was utilized to infer the mutation-scaled effective population size (Θ = 2Nμ, with N = effective population size and μ = mutation rate) and the mutation-scaled effective immigration rate (Μ = m/μ, with m = immigration rate) [61]. Analyses were conducted with a full migration matrix model (Θ and Μ were estimated jointly from the data). The effect number of migrates per generation (Nm) among groups can be calculated by multiplying Θ and Μ together. For the Bayesian approach, a single long chain with slice sampling for the proposal distribution was used according to the recommendation of the author [61]. We performed Migrate with the DNA sequence model. Initially, short runs were performed to estimate Θ and Μ with F ST with a uniform prior for parameters, followed by subsequent runs set using parameters estimated in the short runs. Five independent sets of runs were conducted, each containing one long chain of 5 000 000 steps with a burn-in time of 500 000, a sampling increment of 1000, and an adaptive heating scheme with four chains and temperatures of 1.0, 1.5, 3.0, and 10 000. The results of the independent runs were congruent and the last run was chosen for interpretations.