Genetic Diversity and Population Structure of Acanthochiton rubrolineatus (Polyplacophora) Based on Mitochondrial and Nuclear Gene Markers

: Acanthochiton rubrolineatus (Cryptoplacidae, Neoloricata, Polyplacophora) has a narrow distribution range along the seacoasts of China, the Korean Peninsula and Japan. We collected 238 samples from eight localities along the Chinese coast, and analyzed the genetic diversity and population structure with COI, 16S-rRNA and 28S-rRNA gene sequences. All analyses based on combined sequences of COI and 16S-rRNA suggested that there was evident genetic di ﬀ erentiation between the northern populations (YT, WH, DL, QD, LYG) and southern populations (ZS, YH, XM) of A. rubrolineatus . The haplotype distribution pattern and genetic diversity based on 28S-rRNA sequences also supported the genetic divergence between the two groups. Both groups had experienced population expansion after the ice age of Pleistocene, and an additional population bottleneck had happened in the southern group in recent history, which led to low genetic diversity of mitochondrial DNA and abnormally high diversity of nuclear DNA in this group. Our results suggested that the protection on A. rubrolineatus is necessary, and the northern and southern group should be protected separately.


Introduction
Polyplacophora (Mollusca) contains more than 940 species and 430 fossil species of chitons [1,2], which are oval, dorsoventrally compressed and bilaterally symmetrical marine invertebrates with eight valves in their dorsum [3]. Chitons are a relatively primitive group among mollusks [4], but their shape and living habits have not changed significantly in the past 500 million years [5,6]. Therefore, chitons are called as "living fossils", and they have great value in studies of speciation and the evolution of mollusks [7]. Many investigations about this group have been focused on the origin and phylogenetic relationships among different species [7][8][9][10][11], population structure [12], morphology [13,14] and embryonic development [15][16][17][18]. In addition, as chitons are common species in littoral rocky coasts in their distribution regions, they can serve as good indicators of the local environment. Some species of chitons can be eaten as food and have great medicinal interest [19][20][21][22][23]. The radular teeth of chitons, containing large amounts of magnetite in the form of nanoparticles, were found to be interesting research materials of biomineralization, magnetic materials and natural nanomaterials [24][25][26][27].
Acanthochiton rubrolineatus (Cryptoplacidae, Neoloricata, Polyplacophora) is a species that only distributes along the seacoasts from the Bohai Sea to the East China Sea, and the Korean Peninsula and Japan [3,28]. The three red dark fringes on the intermediate valves, the wide girdle scales, and the granule shape bulges on scales' surface are distinctive morphological features of A. rubrolineatus [3]. Studies on this species have been focused on its great economical values on new magnetic materials Diversity 2020, 12 and natural nanomaterials [24,25] and biopharmaceuticals [20][21][22][23]29]. The concern is that the wild resources of A. rubrolineatus are declining sharply in Chinese coasts due to coast overexploitation and ocean pollution [30]. Our field survey showed that this species has disappeared in some previously recorded localities. Wild resource protection is the only way to conserve A. rubrolineatus, because this species has not been artificially cultured so far [31]. A comprehensive understanding on the genetic diversity and population structure of this species is necessary for the establishment of protection measures. Population genetic diversity and structure can reveal the evolutionary history of the species [32,33]. Mitochondrial gene markers, especially the cytochrome c oxidase subunit I gene (COI) and 16S ribosomal RNA gene (16S-rRNA), are commonly used in genetic diversity analyses of marine invertebrates [34][35][36][37][38]. Now, the genetic diversity and population structure of A. rubrilineatus in the whole distribution region is still unclear, though samples from Bohai Bay have been analyzed with COI sequences [39].
In this study, we collected A. rubrolineatus samples from eight localities along the Chinese coast, which generally represent the whole distribution region of this species in China. We analyzed the population genetic diversity and structure using both mitochondrial (COI, 16S-rRNA) and nuclear (28S-rRNA) gene markers. The results will be greatly helpful for wild resource conservation and fisheries management of A. rubrolineatus.

Sample Collection
A total of 238 individuals of A. rubrolineatus were collected from eight localities across the coasts of China sea (Dalian, DL; Yantai, YT; Wehai, WH; Qingdao, QD; Lianyungang, LYG; Zhoushan, ZS; Yuhuan, YH; Xiamen, XM) (Table 1, Figure 1). All individuals were live trapped and identified based on morphological characters [40]. The abdomen muscle tissues were placed in absolute alcohol and transported to the laboratory, and then were stored at −80 • C until DNA extraction. and Japan [3,28]. The three red dark fringes on the intermediate valves, the wide girdle scales, and the granule shape bulges on scales' surface are distinctive morphological features of A. rubrolineatus [3]. Studies on this species have been focused on its great economical values on new magnetic materials and natural nanomaterials [24,25] and biopharmaceuticals [20][21][22][23]29]. The concern is that the wild resources of A. rubrolineatus are declining sharply in Chinese coasts due to coast overexploitation and ocean pollution [30]. Our field survey showed that this species has disappeared in some previously recorded localities. Wild resource protection is the only way to conserve A. rubrolineatus, because this species has not been artificially cultured so far [31]. A comprehensive understanding on the genetic diversity and population structure of this species is necessary for the establishment of protection measures. Population genetic diversity and structure can reveal the evolutionary history of the species [32,33]. Mitochondrial gene markers, especially the cytochrome c oxidase subunit I gene (COI) and 16S ribosomal RNA gene (16S-rRNA), are commonly used in genetic diversity analyses of marine invertebrates [34][35][36][37][38]. Now, the genetic diversity and population structure of A. rubrilineatus in the whole distribution region is still unclear, though samples from Bohai Bay have been analyzed with COI sequences [39].
In this study, we collected A. rubrolineatus samples from eight localities along the Chinese coast, which generally represent the whole distribution region of this species in China. We analyzed the population genetic diversity and structure using both mitochondrial (COI, 16S-rRNA) and nuclear (28S-rRNA) gene markers. The results will be greatly helpful for wild resource conservation and fisheries management of A. rubrolineatus.

Sample Collection
A total of 238 individuals of A. rubrolineatus were collected from eight localities across the coasts of China sea (Dalian, DL; Yantai, YT; Wehai, WH; Qingdao, QD; Lianyungang, LYG; Zhoushan, ZS; Yuhuan, YH; Xiamen, XM) (Table 1, Figure 1). All individuals were live trapped and identified based on morphological characters [40]. The abdomen muscle tissues were placed in absolute alcohol and transported to the laboratory, and then were stored at −80 °C until DNA extraction.

PCR Amplification and Sequencing
The genomic DNA of each individual was extracted from muscle tissue with a TIANamp Marine Animals DNA Kit (Catalog number: DP324-03, TIANGEN, Beijing, China) according to the manufacturer's instructions, and was adjusted to 50 ng/μL. The gene segments of COI, 16SrRNA and

PCR Amplification and Sequencing
The genomic DNA of each individual was extracted from muscle tissue with a TIANamp Marine Animals DNA Kit (Catalog number: DP324-03, TIANGEN, Beijing, China) according to the manufacturer's instructions, and was adjusted to 50 ng/µL. The gene segments of COI, 16SrRNA and 28SrRNA were amplified respectively with a polymerase chain reaction (PCR) technique using the primers listed in Table S1. The PCR reaction system was same as that described in Gong et al. [37]. The PCR procedures were listed in Table S1. A negative control without genomic DNA was included in each round of PCR to check for contamination, and all negative controls obtained no products. The PCR products were checked on 1.5% agarose gels stained by ethidium bromide and were visualized by ultraviolet transillumination (UVP, Upland, CA, USA). Then, the products were purified with a gel extraction kit (Sangon BioMedical, Shanghai, China) and sequenced (both directions) with an ABI 3730 XL automatic sequencer (Perkin-Elmer, Waltham, MA, USA) using an ABI PRISM BigDye Terminator Cycle Sequencing Ready Reaction kit (with AmpliTaq DNA polymerase FS, Applied Biosystems, Foster City, CA, USA).

Data Analysis
The sequences of COI (652 bp) and 16S-rRNA (533 bp) segments from 238 individuals ( Table 1) and sequences of 28S-rRNA (326 bp) segments from 235 individuals (27 individuals in WH) were obtained separately. All sequences were deposited in GenBank (accession numbers: MT338576 -MT338813 for COI, MT331859 -MT332096 for 16S-rRNA, MT326640 -MT326874 for 28S-rRNA). In the analyses, two mitochondrional genes (COI and 16S-rRNA) were combined. All sequences were aligned using the ClustalX 2.0 [41,42] under default settings with manual corrections as implemented in the software BioEdit v7.0.9 [43], and then were adjusted by visual inspection and BLAST analyses in the NCBI database. All haplotypes were defined by DNASP5.10.01 [44], and the genetic diversity parameters (the number of haplotypes-h, haplotype diversity-Hd, number of segregating sites-S and nucleotide diversity-π), were also calculated using DNASP 5.10.01. To evaluate historical demography, the Fu's Fs statistics and Tajima's D [45] were calculated using the Arlequin 3.5.13 [46] with 1000 permutation to test for neutrality. The values of fixation index (F ST ) between different populations and the analyses of molecular variance components (AMOVA) were also calculated by the Arlequin 3.5.13. Kimura 2-parameter model distances were calculated using Mega 6.0.6 [47]. The expansion times of the southern and northern groups were inferred using a Bayesian skyline plot (BSP) model implemented in BEAST v2.6.2 [48]. The appropriate substitution model of sequences was chosen by Modeltest version 3.7 [49,50]. Coalescent genealogies of BSP were constructed with a lognormal relaxed clock model for a total of 10 million generations by sampling every 1000 generations. The mutation rate of COI sequences is 2%/Mya (the mutation rate of COI in mollusk) [51,52], the appropriate substitution model is TPM3+G, and G is 0.1850. For 16S rRNA sequences, the mutation rate is 1.2%/ Mya [51,52], the appropriate substitution model is GTR +G, and G is 0.6420.
Phylogenetic trees were constructed using Bayesian inference and maximum likelihood methods to explore the relationship among different haplotypes defined from eight geographical localities. The phylogenetic analyses were based on combined sequences of COI and 16S rRNA and sequences of 28S rRNA, respectively. Acanthochitona crinite (COI and 16S-rRNA: MG680027) and Liolophura japonica (28S-rRNA: AY377683.1) were used as outgroup in analyses with two sets of data. Bayesian analyses were performed using MrBayes 3.2.7 [53] with an appropriate substitution model of sequence that was chosen by Modeltest v3.7 [49,50]. The best substitution model for combined sequences of COI and 16S-rRNA (GTR+G) and 28S-rRNA (K80+G) was estimated separately. For the Bayesian procedure, four independent Markov chains were run for 10 million generations by sampling one tree per 1000 generations and allowing adequate time for convergence. The first 2500 trees (25%) were discarded as part of a burn-in procedure (that was determined by checking for the likelihood of being stationary); the remaining 7500 sampling trees were used to construct a 50% majority rule consensus tree. Two independent runs were carried out to provide additional confirmation of the convergence of the Bayesian posterior probabilities (BPP) distribution. The maximum likelihood (ML) tree with 100 bootstrap pseudoreplicates was constructed online with PhyML 3.0 [54] (http: //www.atgc-montpellier.fr/phyml/). The model of evolution chosen by Modeltest v3.7 was applied, but PHYML was allowed to estimate parameter values independently.
A median-joining network was constructed using Network 5.0.1.1 [55] to infer relationships among 92 haplotypes defined by combined sequences of COI and 16S rRNA.

Population Structure and Genetic Diversity Based on Mitochondrial Gene Markers
A total of 92 haplotypes and 249 polymorphic sites (156 singleton variable sites and 93 parsimony informative sites) were defined by 238 combined sequences of COI and 16S rRNA (1185 bp) ( Table 2, Table S2). The Bayesian tree and ML tree, based on sequences of 92 haplotypes, presented consistent topology, which showed clearly that all haplotypes were clustered into two lineages ( Figure 2). All haplotypes from northern China (YT, WH, DL, QD, LYG) gathered into a lineage, and all haplotypes from southern China (ZS, YH, XM) grouped into another lineage. Among 92 haplotypes, 71 haplotypes were private (accounting for 77.2%) (Table S2). Five northern populations shared two haplotypes (H6 and H11), and three southern populations shared one haplotype (H54). In the median joining network, two radiation clades corresponding to the northern and southern group were presented ( Figure 3). In the northern clade, the haplotype H6, shared by five populations, was the dominant haplotype. In the southern clade, the haplotype H54, shared by three populations, was located in the radiation center. The correlations among haplotypes reflected in the network ( Figure 3) were generally consistent with the relationships showed in the phylogenetic tree ( Figure 2).
For all samples, the haplotype diversity (Hd) and nucleotide diversity (π) were 0.889 and 0.03310, respectively. For different populations, the Hd ranged from 0.648 (LYG) to 0.959 (DL), and the π ranged from 0.00152 (ZS) to 0.01059 (QD) ( Table 2). Because two apparent groups were suggested by the phylogenetic tree, the genetic diversity for the two groups were also calculated. For the northern group (DL, YT, WH, QD, LYG), the Hd was 0.858 and the π was 0.00871; for the southern group (ZS, YH, XM), the Hd was 0.639 and the π was 0.00379 (Table 2).    Figure 1, and circle size is proportional to the number of individuals with the same haplotype. The distribution of all haplotypes was also given in Table S2.
For all samples, the haplotype diversity (Hd) and nucleotide diversity (π) were 0.889 and 0.03310, respectively. For different populations, the Hd ranged from 0.648 (LYG) to 0.959 (DL), and the π ranged from 0.00152 (ZS) to 0.01059 (QD) ( Table 2). Because two apparent groups were suggested by the phylogenetic tree, the genetic diversity for the two groups were also calculated. For the northern group (DL, YT, WH, QD, LYG), the Hd was 0.858 and the π was 0.00871; for the southern group (ZS, YH, XM), the Hd was 0.639 and the π was 0.00379 (Table 2).
In the neutrality tests, except the Tajima's D value of DL, the other Tajima's D and Fu's Fs test values for all eight populations were negative ( Table 2). The P values for five populations (QD, LYG, ZS, YH, XM) were smaller than 0.05. When the tests were carried out based on northern and southern groups, all of the Tajima's D and Fu's Fs test values were negative with P values smaller than 0.02.
Based on the combined sequences of mitochondrial gene segments, the FST values among populations from northern region ranged from −0.02071 to 0.13911, and the FST values among populations from southern China ranged from −0.02332 to −0.00791 (Table 3). However, the FST values between northern and southern populations ranged from 0.59467 to 0.74835. The Kimura 2-parameter model distance among populations from the northern region ranged from 0.00604 to 0.01489, and the distance among the populations from southern China ranged from 0.00105 to 0.00579. Between the northern and southern populations, the K2P distance ranged from 0.06106 to 0.07025 (Table 3).
AMOVA analyses based on eight localities showed that 53.19% of the genetic variance of A. rubrolineatus came from different populations (Table S3). When the analysis was carried out based on two groups (the northern and southern group), the results suggested that 66.78% of the genetic variance existed between two groups (Table S3).  Figure 1, and circle size is proportional to the number of individuals with the same haplotype. The distribution of all haplotypes was also given in Table S2. In the neutrality tests, except the Tajima's D value of DL, the other Tajima's D and Fu's Fs test values for all eight populations were negative ( Table 2). The P values for five populations (QD, LYG, ZS, YH, XM) were smaller than 0.05. When the tests were carried out based on northern and southern groups, all of the Tajima's D and Fu's Fs test values were negative with P values smaller than 0.02.
Based on the combined sequences of mitochondrial gene segments, the F ST values among populations from northern region ranged from −0.02071 to 0.13911, and the F ST values among populations from southern China ranged from −0.02332 to −0.00791 (Table 3). However, the F ST values between northern and southern populations ranged from 0.59467 to 0.74835. The Kimura 2-parameter model distance among populations from the northern region ranged from 0.00604 to 0.01489, and the distance among the populations from southern China ranged from 0.00105 to 0.00579. Between the northern and southern populations, the K2P distance ranged from 0.06106 to 0.07025 (Table 3).  (Table S3). When the analysis was carried out based Diversity 2020, 12, 159 7 of 14 on two groups (the northern and southern group), the results suggested that 66.78% of the genetic variance existed between two groups (Table S3).

Population Structure and Genetic Diversity of A. rubrolineatus Inferred from 28S-rRNA Sequences
Only 20 haplotypes and 23 polymorphic sites (10 singleton variable sites and 13 parsimony informative sites) were defined by 235 28S-rRNA sequences (342 bp) from eight localities (Tables S4  and S5). The phylogenetic tree, based on 28S-rRNA sequences (Figure 4), presented different topology from that of the mitochondrial gene markers (Figure 2). There was no distinct lineage corresponding to the northern or southern group, and some haplotypes from the northern and southern regions mixed in one clade (Figure 4).

Population Structure and Genetic Diversity of A. rubrolineatus Inferred from 28S-rRNA Sequences
Only 20 haplotypes and 23 polymorphic sites (10 singleton variable sites and 13 parsimony informative sites) were defined by 235 28S-rRNA sequences (342 bp) from eight localities (Tables S4  and S5). The phylogenetic tree, based on 28S-rRNA sequences (Figure 4), presented different topology from that of the mitochondrial gene markers (Figure 2). There was no distinct lineage corresponding to the northern or southern group, and some haplotypes from the northern and southern regions mixed in one clade (Figure 4). Among 20 haplotypes, eight haplotypes were private (account for 40.0%) (Table S5). Five northern populations shared two haplotypes (S2 and S4), and three southern populations shared six haplotypes (S7, S9, S10, S12-S14). The Hd and π values of all samples were 0.830 and 0.01350, Among 20 haplotypes, eight haplotypes were private (account for 40.0%) (Table S5). Five northern populations shared two haplotypes (S2 and S4), and three southern populations shared six haplotypes (S7, S9, S10, S12-S14). The Hd and π values of all samples were 0.830 and 0.01350, respectively (Table S4). For different populations, the Hd ranged from 0.614 (YT and LYG) to 0.834 (ZS and YH), and the π ranged from 0.00316 (DL) to 0.01975 (ZS). The Hd and π values of the southern populations were greater than those of the northern populations (Table S4).

Significant Genetic Divergence Between Northern and Southern Groups of A. rubrolineatus
Among 92 haplotypes defined by the mitochondrial gene markers, no haplotype was shared by northern and southern populations (Table S2). In the phylogenetic trees, all 58 haplotypes from northern populations (YT, WH, DL, QD, LYG) grouped into one lineage, and 34 haplotypes from  (Figure 2). In the median joining network, haplotypes from northern and southern regions also distributed in two radiation clades, respectively ( Figure 3).
Mitochondrial COI and 16S rRNA sequences were extensively used to assess the genetic diversity of marine invertebrate species [28,37,[56][57][58][59][60][61]. Combining previous reports, the mitochondrial DNA diversity of A. rubrolineatus was relatively high. However, the diversity level of the northern group and southern group was very different, and the overall diversity of the northern group was significantly higher than that of the southern group ( Table 2).
The pairwise F ST value is a popular parameter to estimate gene flow and genetic differentiation between different populations. When the F ST values are greater than 0.05, genetic differentiation exists among populations; if the F ST values are greater than 0.25, the genetic differentiation among populations is significant [62][63][64]. Based on the mitochondrial gene markers, the F ST values among populations from the southern region were negative (−0.02332 to −0.00791) ( Table 3), which indicated that there was no genetic divergence within the southern group. Except DL, the F ST values among the other four northern populations were smaller than 0.05. However, the F ST values between the northern and southern populations were all significantly greater than 0.25 (0.59467 to 0.74835) ( Table 3).
The Kimura 2-parameter model distance between the northern and southern populations (0.06106-0.07025) were significantly greater than those within the northern group (0.00604-0.01489) and southern group (0.00105-0.00579) ( Table 3). The AMOVA analyses also showed that most of the genetic variance (66.78%) existed between two groups (Table S3).
The results of all above analyses, based on combined sequences of mitochondrial gene segments, suggested that gene flow between the northern and southern group of A. rubrolineatus was very limited, and significant genetic divergence had been accumulated between two groups (Figures 2 and 3; Tables 2  and 3; Tables S2 and S3). Our results indicated that the classification of A. rubrolineatus need to be reconsidered, though there was no report about the reproductive isolation between the northern and southern groups of this species.
Because of the low mutation rate of 28S-rRNA, few haplotypes were defined by sequences of this gene (Tables S2 and S5), and the phylogenetic tree based on sequences of this gene did not show two distinct geographic lineages ( Figure 4). However, the haplotype distribution pattern and genetic diversity also presented evidence for obvious divergence between northern and southern groups. Two groups shared no 28S-rRNA haplotype (Table S5), and the genetic diversity of southern populations were greatly higher than those of northern populations (Table S4).
There are several possible reasons to explain the great genetic divergence between the northern and southern groups. First, most adults of chitons live on reefs, and it is difficult for them to migrate a long distance. Moreover, the planktonic larval of chitons have limited capacity for long-distance dispersal [65][66][67], because most chiton species have only lecithotrophic (non-feeding) larvae and their planktonic larval stage is very brief [68]. Second, the Yellow Sea and East China Sea were separated during the subglacial periods of the Pleistocene Ice Age because of the falling water level [69], which would hinder the gene flow between the northern and southern group completely. During the interglacial stage, there are many small and persistent convection and eddy currents near the coasts of China Sea [70], which will limit the planktonic larva of chitons to nearshore areas and furtherly hinder their long-distance migration. DL and the other four northern populations were isolated by the Bohai Sea (Figure 1). The greater F ST values between DL and the other four populations within the northern group (Table 3) can be well explained by the limited migration capacity of the adults and planktonic larvae of chitons. Third, Bohai Bay, the Yellow Sea (monsoon climate of medium latitudes) and East China Sea (subtropical monsoon climate) belong to different climatic types, and their annual temperature change characteristics, the annual precipitation, and the water salinity and pH values are very different [71]. These factors lead to the differences of growth rate, individual size and reproductive period between the northern and southern populations of A. rubrolineatus. Rising temperature and acidification of the seawater would increase feeding amount due to the increase of metabolic demand [72,73]. The body sizes of individuals from southern China are generally larger than those of samples from northern China. In addition, there is a distribution gap of A. rubrolineatus at the southern border of the Yellow Sea (because about 900 km of the coasts in Jiangsu Province is mudflat without rock), which lead to the absence of step-stone populations between the northern and southern group. Last, the different natural climate and coastal ecological environment produced different natural selective pressures on the northern and southern populations, which may furtherly facilitate the accumulation of genetic divergence.
Considering the apparent genetic divergence and different genetic diversity between the two groups, we proposed that the northern group and southern group of A. rubrolineatus should be protected separately.

Historic Demography of A. rubrolineatus
In general, the genetic diversity of A. rubrolineatus is characterized by high haplotype diversity and low nucleotide polymorphism ( Table 2, Table S4). Many marine invertebrate species that have experienced population decline and following expansion showed this type of genetic diversity [37,58,59], because haplotype diversity may increase rapidly in a short period of time by accumulating mutations, while nucleotide polymorphism cannot return to the original level in the same time [74].
A total of 68 individuals in five northern populations shared mitochondrial DNA haplotypes H6 and H11, and 125 individuals shared 28S-rRNA haplotypes S2 and S4 (Tables S2 and S5). There were 53 individuals in three southern populations that shared the mitochondrial haplotype H54, and a total of 80 individuals shared six 28S-rRNA haplotypes (S7, S9, S10, S12-S14) (Tables S2 and S5). The extensively shared haplotypes among different populations is another sign of population expansion.
Tajima's D values and Fu's Fs tests are frequently used parameters in neutrality tests to examine recent population expansion [75]. Negative Tajima's D value with a P value smaller than 0.05 indicates that the sequence contains more nucleotide changes than the neutral evolution model [76]. Based on combined sequences of COI and 16S rRNA segments, the Tajima's D and Fu's Fs test values for five populations (QD, LYG, ZS, YH, XM) were negative with P values smaller than 0.05, and the values were negative with P values smaller than 0.02 when the estimation was carried out for two (northern and southern) groups ( Table 2). The patterns of genetic diversity and haplotype distribution, as well as the results of the neutrality test, suggested that both groups of A. rubrolineatus had experienced population expansions.
We calculated the expansion times for two groups using the extended Bayesian skyline plot [48] based on COI and 16S rRNA sequences, respectively. Based on COI sequences, the expansion time was about 0.349 Mya and 0.205 Mya for the northern group and southern group, respectively. Based on 16S rRNA sequences, the expansion time was about 0.280 Mya and 0.160 Mya for the northern group and southern group, respectively. The estimated time was corresponding to the interglacial stage after the Dagu-Lushan glaciation of middle Pleistocene [77].
In the quaternary glacial epoch, there were at least four subglacial and interglacial periods in China. Previous studies on population structures of many marine invertebrate species from Chinese seacoasts [37,[78][79][80][81] showed that Pleistocene glaciation was the main driving force to shape the genetic structures of marine mollusks in China. The dramatical shrinkage of the China Sea and temperature decrease in the subglacial periods led to the population decline or extinction of many marine invertebrate species in some regions. During the last glacial maximum (21,000 BP-15,000 BP), Bohai Bay and the Yellow Sea almost disappeared [69], the East China Sea became the shelter of most marine mollusk species. As a primitive mollusk species that has existed for hundreds of million years, A. rubrolineatus must have experienced several times population decline and following population expansion during the Pleistocene Ice Age. Under this condition, the expansion time estimated from sequences of a few genes cannot reflect the detailed population history of A. rubrolineatus, and the error of calculation is easy to understand.
Previous studies on population structures of other marine mollusk species (such as Atrina pectinate, Coelomactra antiquata, Eriocheir sensu stricto, Urechis unicinctus) did not find obvious genetic divergence between northern and southern populations [37,[78][79][80]. Therefore, it was proposed that the existing populations of many marine mollusk species in Bohai Bay and the Yellow Sea migrated from the East China Sea by population expansion after the last glacial maximum [81]. However, the situation in A. rubrolineatus was very different, the great genetic divergence between northern and southern group and the different expansion time estimated by mitochondrial genes indicated that the existing populations of A. rubrolineatus along Bohai Bay and the Yellow Sea are descendants of native survivors in the ice age.
We should notice that the genetic diversity of the two groups inferred from mitochondrial and nuclear gene segments were very different ( Table 2, Table S4). Considering the lower mutation rate of 28S-rRNA than mitochondrial genes, the lower genetic diversity of 28S-rRNA than that of mitochondrial genes in northern groups should be a normal situation. The low genetic diversity of mitochondrial DNA and abnormally high diversity of nuclear DNA in the southern group indicated that an additional population bottleneck had happed in recent history. Strong typhoons and storm surges often attack coasts of the East China Sea, while the coasts of Bohai Bay and the Yellow Sea are rarely attacked by such disasters [82]. Strong typhoons and storm surges always destroy coastal rocks and reefs severely, and habitat destruction will lead to the population decline of A. rubrolineatus in southern China. The bottleneck effect, combining the genetic drift effect, result in the special genetic diversity pattern in southern populations.

Conclusions
A. rubrolineatus is a species of chitons with a narrow distribution area. This species has great scientific and economical values, but its wild resources are declining sharply. A comprehensive understanding of the genetic diversity and population structure of this species is necessary for the establishment of protection measures. In the present study, 238 individuals of A. rubrolineatus were collected from eight localities along Chinese seacoasts, and the population genetic diversity and structures were analyzed with both mitochondrial (COI, 16S-rRNA) and nuclear (28S-rRNA) gene markers. The results showed that there was great genetic differentiation between the northern and southern group of A. rubrolineatus and both groups had experienced population expansion. An additional population bottleneck had happened in the southern group in recent history, which led to the different patterns of genetic diversity in the two groups. We proposed that the northern group and southern group of A. rubrolineatus should be protected separately.