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

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.

Table 2. Sampling sites and summary of diversity indices of Siphonaria japonica. For each locality, individual numbers (n), haplotype diversity (h) and nucleotide diversity (?) for both mitochondrial and nuclear DNA are listed

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. toreuma20].

Table 3. Pairwise genetic distances (?_ST) among locations of Siphonaria japonica. Pairwise ?_ST values for nuclear sequence ITS and mitochondrial sequence COI of Siphonaria japonica are given in the upper and lower diagonals, respectively. *P??0.05; **P??0.01, ***P??0.001

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).

Fig. 2. Maps of gene flow and median-joining network (Insert) for each studied species. Map
of the China coast showing sampling sites and estimated number of migrants per generation
between groups (Nm) for each species. All sampling sites of each species are divided into Yellow Sea
(YS) group and East plus South China Seas (ESCS) group according to the Yangtze River
Estuary. The black dots represent sampling localities for each species. The blue circle
and red circle represent the relative magnitude of ?_YS (YS group) and ?_ESCS (ESCS group) respectively. The blue and red arrows represent the direction and relative
magnitude of gene flow estimates between groups. Inset: Haplotype network obtained
with mitochondrial DNA for each species and nuclear DNA for S. japonica. Each circle represents a single haplotype and sizes are proportional to the number
of individuals possessing these haplotypes. Haplotypes of each species are divided
into the Yellow Sea (YS) group and the East and South China Seas (ESCS) group separated
by the Yangtze River estuary. Colors represent group where haplotypes were detected
(blue, YS group; red, ESCS group). The small white circles indicate hypothetical missing
haplotypes. Reference: Cellana toreuma, Dong et al. 20]; Sargassum horneri, Hu et al. 36]; Atrina pectinata, Liu et al. 37]

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 expansion (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.

Fig. 3. Mismatch distributions of the observed haplotypes for each species. The histograms
are the observed frequencies of pairwise divergences and the line refers to the expectation
under the sudden population expansion model. All COI sequences of S. japonica were divided into two putative 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]

Fig. 4. Bayesian skyline plots (BSP) for Siphonaria japonica. All COI sequences of S. japonica were divided into two groups: northern group (a) and southern group (b), according to the phylogenetic analyses (See the results in the phylogenetic analysis).
The black line represents median population estimates; the upper and lower limits
of light blue shading represent the 95 % confidence intervals. The gray rectangles
represent the glacial stages of various oxygen isotope stages 42], 43] in (c)

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
were unimodal (Fig. 3c, d, e) and the time of demographic expansion for C. toreuma, S. horneri and A. pectinata were 502.5 ka (90 % CI: 70.7-538.2 ka) BP, 119.3 ka (90 % CI: 47.7-195.5 ka) BP and
170.1 ka (90 % CI: 133.7-421.7 ka) BP respectively (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 rocky-shore 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 the mitochondrial sequences, there is southward gene flow from YS
to ESCS (Nm_(YS-ESCS)?=?5.79-24.93), but little gene flow northwards from ESCS to YS (Nm_(ESCS-YS)?=?0.00-8.14). The ITS analyses, however, indicated there was almost no gene flow in a southward direction
from YS to ESCS (Nm_(YS-ESCS)?=?0.00-4.55), compared with a high gene flow from ESCS to YS (Nm_(ESCS-YS)?=?18.37-138.29).

Cellana toreuma, Sargassum horneri and Atrina pectinata

Re-analysis of published data showed that significant higher southward gene flow could
be observed in C. toreuma (Nm_(YS-ESCS)?=?43.41-359.13, Fig. 2c) and A. pectinata (Nm_(YS-ESCS)?=?34.87-3620.72, Fig. 2e). On the contrary, there is almost no gene flow in a southward direction (Nm_(YS-ESCS)?=?0.00-31.57, Fig. 2d) in S. horneri, while large gene flow can be observed from ESCS to YS (Nm_(ESCS-YS)?=?59.93-3629.67).