Phylogeography and Taxonomic Revision of the Pen Shell Atrina pectinata Species Complex in the South China Sea

Phylogeographic studies contribute to addressing questions regarding the geographic patterns and evolutionary scenarios within and among species and also shed light on the taxonomic status of widely distributed species complexes. The pen shell Atrina pectinata species complex is a widely distributed and economically important bivalve in the northwestern Pacific. Previous phylogeographic studies have identified four genetically distinct cryptic species within the A. pectinata species complex along the coast of China, of which three cryptic species were distributed in the South China Sea. However, less attention has been given to their identification and delimitation. Herein, we report the phylogeography and taxonomic revision of the A. pectinata species complex in the South China Sea using DNA sequence data from mitochondrial cytochrome oxidase I (COI) and 16S ribosomal RNA (16S rRNA), and nuclear 28S ribosomal RNA (28S rRNA) gene markers. Using a combination of phylogenetic and DNA-based species delimitation analysis methods, we found strong support for four genetically valid species in the A. pectinata species complex and defined them as A. japonica, A. lischkeana, Atrina sp., and A. pectinata based on our results as well as on previous morphological and genetic studies. A demographic historical analysis showed that all three species in the South China Sea had populations that were relatively stable over time and then subjected to sudden expansion during the late Pleistocene (60,000–90,000 years ago). These results provide new insights into the systematics and evolution of the A. pectinata species complex and have important conservation and management implications.


INTRODUCTION
Species identification and delimitation are fundamentally important within the fields of biology, biogeography, ecology, and conservation, as species are the fundamental units for biodiversity quantification and management (Agapow et al., 2004;Bickford et al., 2007;Liu et al., 2020). Historically, species identification and delimitation have been performed by taxonomists using morphological characteristics; these processes require a high level of expertise and are time consuming. However, many diverse groups are morphologically cryptic and contain many undescribed taxa, and the existing taxonomic literature is conflicting (Bickford et al., 2007). Molecular analytical approaches (e.g., DNA barcoding) have been developed to identify and delimit species and are now used to resolve many taxonomic problems associated with different animals and plant species (Wiens and Servedio, 2000;Goldstein and DeSalle, 2011;Previšić et al., 2016;Sun et al., 2016;Liu et al., 2020). Several methods have been provisionally designed to assess species status based on DNA sequences, including automatic barcode gap discovery (ABGD) (Puillandre et al., 2012), the generalized mixed Yule coalescent model (GMYC) (Fujisawa and Barraclough, 2013),Bayesian Poisson tree processes (bPTP) (Zhang et al., 2013), and Bayesian phylogenetics and phylogeography (BPP) (Yang and Rannala, 2010;Yang, 2015).
The family Pinnidae Leach, 1819, includes approximately 50 species of commercially important subtidal and coastal species that inhabit tropical and temperate coastal oceanic waters worldwide (Huber, 2010;Lemer et al., 2014). The taxonomy of the family has been subjected to a number of revisions since the early 20th century, and the overall number of valid species ranges from 21 to 175 in various reports (Lemer et al., 2014). Atrina pectinata Linnaeus, 1767, one of the most diverse Pinnidae species, is a large wedge-shaped marine bivalve distributed along the coasts of Indo-Pacific countries. The taxonomy of A. pectinata in East Asia has been complicated because of its shell morphology and lack of unified morphological characteristics (Hashimoto et al., 2018). Since the first description of this species in 1767, the nomenclature of A. pectinata has been confusing: Rosewater (1961) assigned 12 distinct Indo-Pacific Atrina species to a single A. pectinata species, Wang (1997) assigned 11 synonymous names to A. pectinata, and Huber (2010) assigned only three synonymous names to A. pectinata.
In China, five morphologic A. pectinata forms exist within the current concept of A. pectinata (hereafter the "A. pectinata species complex") (Wang, 1997;Yu et al., 2000;Liu et al., 2011). However, debate as to whether the present species should be divided to two or more species or subspecies have lasted for years (Yu et al., 2004;Huber, 2010;Liu et al., 2011;Lemer et al., 2014). Liu et al. (2011) identified five genetically distinct cryptic species (lineage 1-lineage 5) within the A. pectinata species complex along the coast of China. Lineage 1 (nonscaly form) is mainly distributed along the coasts of North China, South Korea and Japan and has been identified as A. japonica (Huber, 2010) and A. pectinata japonica (Kim et al., 2017). Lineage 2-lineage 5 are often sympatric in southern Chinese coastal areas and correspond to four morphological types (thorny, green, yellow, and scabrous pen shells) . Analyses of their morphological characteristics, isoenzymes, random amplified polymorphic DNA (RAPD) and mitochondrial COI sequences have provided evidence of obvious differences among these forms in the South China Sea Yu et al., 2000Yu et al., , 2004Yu et al., , 2019Liu et al., 2011). This evidence strongly implied that extensive morphological work and molecular work are necessary to clarify the taxonomic status of the A. pectinata species complex and suggest that there is less phenotypic plasticity in A. pectinata than previously thought. However, less attention has been given to their identification and delimitation. Lineage 2 (the thorny form) has been called Ken-type in Japan and is thought to be A. lischkeana (Huber, 2010) or A. pectinata lischkeana (Kim et al., 2017;Hashimoto et al., 2018). The morphological and genetic characteristics of lineage 4 (the yellow form) match those of A. pectinata Linnaeus, which is restricted to Indonesia in the Indian Ocean (Huber, 2010;Lemer et al., 2016). Lineage 5 (the scabrous form) has been identified as A. chinensis through the investigation of morphological and genetic data (Xue et al., 2012). However, lineage 3 (the green form) was only reported as one morphological type of A. pectinata distributed in South China Sea (Yu et al., , 2004. As one of the most ecologically and economically valuable marine bivalves, unstable nomenclature could cause management aberrations, so there is an urgent need for species identification and delimitation to support the protection and management of this important bivalve species. In the present study, we collected 211 individuals of the A. pectinata species complex from the South China Sea and used phylogenetic analyses based on mitochondrial cytochrome oxidase I (mtCOI) and 16S ribosomal RNA (mt16S rRNA), and nuclear 28S ribosomal RNA (nr28S rRNA) gene markers from these samples to answer three major questions: (1) Should these four cryptic species (lineage 1 -lineage 4) within the A. pectinata species complex be recognized as distinct taxonomic entities?
(2) What are the evolutionary relationships among these three cryptic species (lineage 2 -lineage 4) in the South China Sea?
(3) Do the three sympatric cryptic species in the South China Sea (lineage 2 -lineage 4) have similar phylogeographic patterns and demographic histories? Their resolution would provide useful information for definite identification and delimitation of the A. pectinata species complex, and developing appropriate conservation and management strategies for the A. pectinata species complex in China, and help reveal the biogeographic history of the studied region.

Sample Collection and Sequencing
A total of 211 individuals of the A. pectinata species complex were collected from 15 sites in southern China from 2010 to 2019 (Figure 1 and Table 1). Hereafter, we defined lineage 1-lineage 4 described in Liu et al. (2011) as four putative species: A. japonica, A. lischkeana, Atrina sp., and A. pectinata, respectively (Figure 2). Adductor muscle samples were taken from each individual and preserved in absolute ethyl alcohol. Genomic DNA was isolated from the ethanol-fixed adductor muscle tissue using the TIANamp Marine Animals DNA Kit (Tiangen; Beijing, China). PCR amplification of three target gene fragments (nr28S rRNA, mtCOI, and 16S rRNA) was performed using the universal nr28S primers D1F and D6R (Park and Foighil, 2000), mtCOI primers LCO1490 and HCO2198 (Folmer et al., 1994), and mt16S primers 16Sar and 16Sbr (Kessing et al., 1989). The mtCOI gene was amplified in 211 individuals, while the nr28S rRNA segment and mt16S rRNA segment were amplified in a subset of individuals FIGURE 1 | Map of the sampling localities and distribution of the three putative species in the A. pectinata species complex along the coasts of southern China. Sampling sites are labeled with abbreviations in Table 1. The accompanying pie charts with different colors indicate their topological placement among the three putative species for each locality. Dotted circles of sampling localities showed that the sequences of these populations were retrieved from GenBank. Abb., Abbreviation of populations, date of collection, N, number of specimens; n, number of specimens sequenced; Nh, number of haplotypes; and np, private haplotypes, are shown for each population. *Sequences of population were retrieved from GenBank (Liu et al., 2011).

Sequence Analyses
Initial alignments were performed using CLUSTAL X2 (Larkin et al., 2007). The sequences were trimmed to the same length with other published sequences after alignment. Molecular diversity indices such as the haplotype diversity (h) and nucleotide diversity (π) of the nr28S rRNA, mtCOI, and mt16S rRNA sequences were calculated in DnaSP 6.12 (Rozas et al., 2017).

Species Delimitation
We applied the ABGD, bPTP, and GMYC methods to test species delimitations using the mtCOI haplotype dataset of A. pectinata species complex. The ABGD method clusters sequences into hypothetical species based on differences between intraand interspecific distance variations (Puillandre et al., 2012). The ABGD analysis was conducted using the online ABGD server 1 under the default parameters. The input tree required for bPTP analyses was generated using MrBayes 3.2 (Ronquist et al., 2012). The bPTP analysis is implemented in an online web server 2 using default model. For the GMYC model (Bouckaert et al., 2019), the ultrametric tree was constructed in BEAST v2.6.3 (Bouckaert et al., 2019) using a Yule pure birth model tree prior. An uncorrelated relaxed lognormal clock model was used with lognormal relaxed distribution. The MCMC length was 100,000,000 generations with log parameters every 10000 generations. The evaluation of ESS values (>200) and trace files of runs were performed in Tracer v1.6. A maximum clade credibility consensus tree was performed obtained in TreeAnnotator v2.6.3. Single-threshold GMYC analyses were carried out in R studio using the splits packages (Ezard et al., 2009). We used the Bayesian Phylogenetics and Phylogeography (BP&P v3.4) program (Yang, 2015) with the three DNA sequence datasets to validate the number of putative species based on the ABGD, bPTP, and GMYC results. The A10 analysis (species delimitation using a guide tree) was repeated twice, with 1,000,000 rjMCMC generations, a burn-in of 8,000, and sampling intervals of 2.
The guide tree used for species delimitation was obtained from phylogenetic analyses (of NJ, ML, and BI gene trees) (Yang, 2015). The species validation was based on the posterior probability value of the species assignment, which indicated the statistical validity of the species.

Phylogenetic Analyses
The nucleotide sequences obtained in this study and those of other pen shell species used for phylogenetic analyses are available from GenBank, and their accession numbers are shown in Supplementary Table 1. P. bicolor (Pinnidae) was used as the outgroup. The best substitution models for the three gene datasets to be used in the phylogenetic analyses were selected based on the Akaike Information Criterion (AIC) with the program Modeltest 3.7 (Posada and Crandall, 1998). Bayesian inference (BI) trees were estimated using MrBayes 3.2 (Ronquist et al., 2012) for 10,000,000 generations with a sample frequency of 1000 generations and an initial burn-in of 100,000. Maximum likelihood (ML) analysis was performed using the PhyML 3.0 program (Guindon et al., 2010). Neighbor-joining (NJ) analyses were conducted using PAUP 4.0b1.0 (Swofford, 2002). The reliability of the internal branches of the NJ and ML trees was estimated with bootstrap values of 1,000 replicates. In addition, the genetic relationships among the haplotypes of the mtCOI genes were further assessed using a minimum spanning tree (Bandelt et al., 1999) constructed by PopART-1.7 3 . The pairwise sequence divergences and the net average genetic distances among species were calculated with MEGA 6.06 (Tamura et al., 2013) using Kimura's two-parameter model (K2P) (Kimura, 1980). In the absence of a clear fossil or geological record, we used mtCOI divergence rates estimated for molluscan species from 0.7% to 2.4% per million years (Myr) (Liu et al., 2011) to assess conservatively a range of dates for key nodes in the mtCOI phylogenetic tree COI.

Population Genetic Structure
To investigate genetic differentiation among populations within the three putative species in South China Sea, analysis of molecular variance (AMOVA) was performed for the mtCOI data using Arlequin3.5 (Excoffier and Lischer, 2010) with 10,000 permutations. AMOVA analyses were conducted with three groups for Atrina sp., the Guangdong group, Hainan group, and Beibu Gulf group, and with two groups for A. lischkeana and 3 http://popart.otago.ac.nz A. pectinata, the Guangdong group and Beibu Gulf group. The pairwise genetic divergence values between populations of the three putative species in South China Sea were estimated using F ST values for mtCOI sequences with Arlequin3.5 (Excoffier and Lischer, 2010), and the significance was adjusted using sequential Bonferroni corrections (Narum, 2006).

Demographic History
Inferences on the historical demographic history of the three putative species in South China Sea were obtained using neutrality tests, mismatch distributions, and Bayesian skyline plots based on the mtCOI sequences. For the neutrality test, Tajima's D test (Tajima, 1989) and Fu's Fs test (Fu, 1997) were calculated using Arlequin3.5 (Excoffier and Lischer, 2010) with 10,000 permutations. A mismatch distribution was constructed for each geographic population to test the exponential population growth model (Rogers and Harpending, 1992). A goodness of fit test was performed to test the validity of the sudden expansion model using a parametric bootstrap approach based on the sum of square deviations (SSD) between the observed and expected mismatch distributions. The raggedness index (RI), which measures the smoothness of the mismatch distribution, was calculated for each distribution. The demographic expansion parameter (τ) was calculated with Arlequin3.5 (Excoffier and Lischer, 2010). The changes in the effective population size of the three lineages with time were also inferred using the Bayesian Skyline method (Ho and Shapiro, 2011) implemented in the BEAST1.7.5 program (Drummond et al., 2012) with 20 groups. Chains were run for 100 million steps to yield effective sample sizes (ESSs) of at least 200, and the first 10% were discarded as "burn-in" under the HKY substitution model (A. lischkeana and A. pectinata) or the GTR substitution model (Atrina sp.) from Modeltest 3.7 (Posada and Crandall, 1998), a strict molecular clock and a stepwise skyline model. All operators were automatically optimized. The results of the analyses were visualized using Tracer 1.5 (Rambaut and Drummond, 2007).

Species Delimitation
The ABGD method clustered the sequences into four groups that were congruent with the four clades inferred in the phylogenetic analyses, representing four putative species: A. japonica, A. lischkeana, Atrina sp., and A. pectinata. These four groups were inferred for a prior maximal distance P ≥ 0.0599 independent of the relative width of the barcoding gap (X-value) or prior specified P range (Supplementary Figure 1). The bPTP model suggested that the estimated number of species is four (Supplementary Figure 2). Using the GMYC approach, the null model of a single species was rejected for mtCOI genes (likelihood ratios for single threshold model: 63.09 (P = 1.99 × 10 −14 ), and indicated four clusters (Supplementary Figure 3). The results of the Bayesian species delimitation based on the three DNA sequence datasets strongly supported that the four major clades defined by the phylogenetic analyses represented four different species, which was congruent with the ABGD, bPTP, and GMYC analyses. Similar results were obtained using different prior distribution sets for θ and τ 0. In summary, the results of ABGD, bPTP, GMYC, and BPP analyses supported the validity of A. japonica, A. lischkeana, Atrina sp., and A. pectinata as their own taxonomic units instead of as synonyms of A. pectinata.

Sequence Variations
Three haplotypes defined by six polymorphic sites were detected among the 44 nr28S rRNA sequences. Each of A. lischkeana, Atrina sp., and A. pectinata samples had only one 1060-bp-long haplotype ( Table 2). The nr28S rRNA sequences of A. lischkeana were exactly the same as those of A. japonica. However, four variable sites were detected between A. lischkeana and Atrina sp., three variable sites were detected between A. lischkeana and A. pectinata, and six variable sites were detected between Atrina sp. and A. pectinata. The overall haplotype diversity (h) and nucleotide diversity (π) values of the 44 sequences were 0.729 and 0.0036, respectively. A 625-bp fragment of the mtCOI gene sequenced for the 211 pen shells along with 117 sequences obtained from GenBank generated 77 haplotypes ( Table 2). Twenty-three haplotypes were observed for 73 individuals of A. lischkeana, 20 haplotypes were observed for 85 individuals of A. pectinata, and 35 haplotypes were observed for 170 individuals of Atrina sp. Nucleotide variations were found in 131 base loci, of which 107 were parsimony informative in the mtCOI dataset. Totals of 56, 36, and 39 variable sites were detected between the mtCOI sequences of A. lischkeana and those of A. pectinata, Atrina sp. and A. japonica, respectively. Forty-six, 51, and 53 variable sites were detected between Atrina sp. and A. pectinata, between Atrina sp. and A. japonica and between A. pectinata and A. japonica, respectively. The overall h and π values were 0.813 and 0.064, respectively.
The sequences of the mt16S rRNA region amplicons from the 52 pen shell individuals were 458 bp in length. Of the seven haplotypes obtained, A. lischkeana had 3 haplotypes, A. pectinata had 3 haplotypes, and Atrina sp. had 1 haplotype ( Table 2). In the mt16S rRNA dataset, 38 sites were variable, of which 37 sites were parsimony informative. Totals of 24, 19, and 10 variable sites were detected between the mt16S rRNA sequences of A. lischkeana and those of A. pectinata, Atrina sp. and A. japonica, respectively. Twenty-seven, 18, and 24 variable sites were detected between Atrina sp. and A. pectinata, between Atrina sp. and A. japonica and between A. pectinata and A. japonica, respectively. The overall h and π values were 0.729 and 0.036, respectively.

Phylogenetic Analysis
The program Modeltest3.7 selected TRN + I (0.4076), HKY + I (0.5503) +G (0.4487), GTR + I (0.5370) +G as the best models based on the AIC for the nr28S rRNA, mtCOI and 16S rRNA datasets, respectively. Phylogenetic analyses of the three DNA markers conducted using the BI, NJ, and ML procedures produced almost identical results (Figures 3-5). A. lischkeana and A. japonica formed one robust subclade and then formed a robust clade with A. pectinata and Atrina sp., which is a sister clade to other Atrina species. Table 3 presents the K2P genetic distances obtained based on the three gene fragments for the above taxa. The pairwise genetic distances among A. lischkeana, Atrina sp., A. pectinata and A. japonica based on 28S sequences were ranged from   (Liu et al., 2011), divergence events among these genealogical clades occurred approximately 3.4-18.1 million years ago during the Neogene period.

Population Demography
The results of Tajima's D and Fu's Fs tests for A. lischkeana, Atrina sp., A. pectinata based on the mtCOI datasets are shown in Table 2. The output values of both tests were significantly negative for all three species, indicating a possible historical population expansion. Furthermore, mismatch distributions for the three species were unimodal (Figure 6), indicating that each of them has experienced a demographic expansion. The low and non-significant SSD (expect A. pectinata) and RI value coincided with the null hypothesis of sudden expansion model ( Table 2). The Bayesian skyline plots demonstrated demographic expansion in each of the three species (Figure 7). Based on the mtCOI divergence rate of 2.4%/Myr, the population expansion times of the A. pectinata and Atrina sp. linkages occurred approximately 60,000 years ago, while that of A. lischkeana was older than those of the other two lineages, occurring approximately 90,000 years ago.

Population Structure
For A. lischkeana, Atrina sp., A. pectinata, the topology of the minimum spanning tree of haplotypes revealed three lineages that corresponded to the three species and were characterized by star shapes (Figure 8). The haplotype network for each  (Figure 6). For A. lischkeana, the global AMOVA analysis showed that there was a shallow but significant amount of variance at only one  level (among populations within groups, F SC = 0.045, P = 0.049). For A. pectinata, there were significant genetic subdivisions among populations within groups (F SC = 0.063, P = 0.030) and within populations (F ST = 0.061, P = 0.006). For Atrina sp., the global AMOVA analysis indicated no significant genetic differentiation at any level (among groups, among populations within groups, or within populations) ( Table 4).
For A. pectinata, the pairwise F ST values ranged from 0.0476 (the ZJ population and BH population) to 0.0921 (the ZJ population and CT population), and two (both involving the ZJ samples) out of three pairwise F ST values were significant after sequential Bonferroni corrections (Table 5), indicating significant genetic differentiation between Zhanjiang and Beibu Gulf (the CT population and BH population). However, no significant pairwise F ST values were detected after sequential Bonferroni corrections for A. lischkeana or Atrina sp.

DISCUSSION
To reach a conclusion on the taxonomy of the A. pectinata species complex, we analyzed the phylogenetic relationships and species delimitations of the A. pectinata species complex in the South China Sea based on the mtCOI, 16S rRNA, and nr28S rRNA gene fragments. In phylogenetic analyses, specimens of the A. pectinata species complex from the South China Sea formed three clusters, and species delimitation analyses indicated that the three lineages revealed in phylogenetic analyses represented three distinct species. These results strongly implied that there is less phenotypic plasticity in A. pectinata than was previously thought. Finally, our data and analyses confirmed previous inferences (Liu et al., 2011) that the four cryptic species of the A. pectinata species complex are genetically valid, and we defined lineage 1-lineage 4 described in Liu et al. (2011) as A. japonica, A. lischkeana, Atrina sp. and A. pectinata based on our results as well as on previous morphological and genetic studies.

Identification and Distribution of the Four Atrina Species in China
Rosewater (1961)'s masterpiece was the reduction of a dozen distinct Indo-Pacific Atrina species into a single species: A. pectinata. However, Winckworth (1936) and Huber (2010) suggested that A. pectinata is restricted to the Indian Ocean and Indonesia and that A. pectinata Linnaeus was distributed in India. A. pectinata specimens collected from the South China Sea shared common haplotypes of mtCOI, 16S rRNA, and nr28S rRNA gene fragments with A. pectinata specimens collected from the Philippine archipelago (Lemer et al., 2016), indicating that A. pectinata is also distributed in the South China Sea. A. pectinata is tan to brownish, rather fragile, narrow, and weakly ribbed, with occasional small spines on the ventral shell surface (Figure 2C).
A. japonica is a large, dark, olive-green to brownish, rather smooth, weakly ribbed, rather inflated species with a very large, central, round muscle scar (Figure 2A) and is mainly distributed from the Yellow Sea to Japan. It was first described by Reeve (1858) in Japan and then synonymized with A. pectinata by some taxonomists (Winckworth, 1936;Habe, 1953;Rosewater, 1961; Wang, 1997). However, Huber (2010) kept it as a distinct species, and Fuxue and Zhang (1959) kept it as a subspecies of A. pectinata. Our results supported A. japonica as a valid species. The shell of A. lischkeana is rather small, fragile and brownish with radiating scaly ribs ( Figure 2B). We differentiated A. lischkeana and A. japonica based on the presence/absence of scales on the dorsal and posterior shell surfaces. A. lischkeana is mainly distributed in Ariake Bay, the Seto Inland Sea in Japan, South Korea and the South Sea in China. A. lischkeana, as described by Clessin (1861) in Yokohama, was synonymized with A. pectinata by Winckworth (1936); Habe (1953), Rosewater (1961), andWang (1997). However, Huber (2010) kept it as a distinct species, and Fuxue and Zhang (1959) kept it as a subspecies of A. pectinata. Our results supported A. lischkeana as a valid species.
Atrina sp. corresponds to the pen shells with green morphological forms  and with lineage 3 in Liu et al. (2011). Atrina sp. comprises rather small, thin, wedgeshaped smooth and green species (Figure 2D). A BLAST search of the three gene fragments of Atrina sp. specimens in this study in GenBank only obtained similar sequences with the highest specimen matches (99-100%) in the South China Sea, indicating that Atrina sp. may be mainly distributed in the South China Sea.
A. japonica is widely distributed over the central to northern seas in China and is more northerly distributed than A. lischkeana, Atrina sp. and A. pectinata, which are mainly distributed in the South China Sea. In areas such as the Leizhou Peninsula and Beibu Gulf, collections of pen shells from Zhanjiang, Caotan, and Beihai have shown that A. lischkeana, Atrina sp. and A. pectinata are sympatric. In collections from Hainan Island, Atrina sp. is most common species found, with only one recorded A. pectinata individual. Two scenarios may potentially account for this distribution pattern: (1) sympatric speciation, in which the speciation process was unrelated to geographic isolation, or (2) pen shells may have originated in tropical regions and expanded their ranges from southern to northern regions after glaciation, giving rise to new lineages/species (Kim et al., 2017). Recently, sympatric speciation has come to be understood as a major generator of marine biodiversity, but it remains among the most controversial evolutionary processes (Quesada et al., 2007;Bowen et al., 2013;Richards et al., 2019). Nearly all existing case studies of sympatric speciation involve some form of automatic magic traits, such as assortative mating by habitat, along a depth gradient, or environment-induced phenology shifts (Richards et al., 2019). However, there is little information on any magic traits associated with the three sympatric pen shell species considered herein, and genomic, physiological, environmental, morphological and reproductive datasets must be combined to reveal speciation mechanisms. According to the lineage-level divergence rates utilized in the present study, the three species separated from one another approximately 3.4-18.1 million years ago during the Neogene period.

Historical Demography of the Four Atrina Species in China
The late Pleistocene period was punctuated by a series of large glacial-interglacial changes that have been proven to affect the phylogeographic distribution and evolutionary processes of marine species (Maggs et al., 2008;Ni et al., 2014). When ice caps formed during Pleistocene glacial periods, coastal marine organisms either went extinct or survived in isolated refugial habitats characterized by the survival of ancestral relict species (Maggs et al., 2008;Hu et al., 2018). The isolation of species within and among glacial refugia and the timing and mode of expansion when glaciers retreated have been widely recognized as important topics for understanding recent speciation and extinction events (Hu et al., 2018). During Pleistocene glacial periods, the sea level in the East China Sea, the Bohai Gulf, and the Yellow Sea (hereafter considered a single sea with no geographical barrier) declined by approximately 150 m, and the East China Sea shelf was completely exposed as the coastline migrated approximately 1,200 km seaward and then consequently receded into an elongated enclosed sea (with the Okinawa Trough) with an area <1/3 of its present size. Meanwhile, the sea level in the South China Sea became approximately 100∼120 m lower than the present day level, and the sea became a semiclosed sea connected to the Pacific Ocean mainly through the Bashi Strait (Ni et al., 2014). The East China Sea and the South China Sea have been proposed as independent marine refugia where taxa adapted to environmental changes by retaining their ancestral polymorphisms during glaciation events, and the dramatic population expansion times of various marine species in these marginal seas were estimated to occur in a period 120,000∼140,000 years ago assuming an mtCOI divergence rate of 2%/Myr (Ni et al., 2014). The mtCOI was utilized as the stable reference marker in historical demography analyses (Marko et al., 2010;Ni et al., 2014). There was neither an accurate mutation rate nor a clear fossil record available for the species. In the former molluscan studies on demographic history, the mtCOI divergence rates estimated ranged from 0.7 to 2.4% per million years (MY) [e.g., Crassostrea gigas: 2∼2.4%/MY ; Thais clavigera: 0.79%/MY (Guo et al., 2015); Cyclina sinensis: 0.7∼2.4%/MY (Ni et al., 2012); Atrina pectinata: 0.7∼2.4%/MY (Liu et al., 2011); Tegula brunnea: 2.4%/MY (Hellberg and Vacquier, 1999)]. Conversion of population genetic parameters with the relatively fast 2.4% rate was expected to yield conservative demographic results with respect to our conclusions: slower mutation rates would increase divergence time estimates and make inferred population expansion times older (Marko et al., 2010). The neutrality test, mismatch distribution and network analyses results combined with Bayesian skyline plots consistently showed that all three studied species in the South China Sea had populations that were relatively stable over time and then subjected to sudden expansion during the late Pleistocene (60,000-90,000 years ago), assuming an mtCOI divergence rate of 2.4%/Myr. This divergence rate is just a rough approximation, causing our estimates of absolute time and population size values to be unreliable. However, the fold-change in Ne and the relative placement of expansions and declines along the temporal axis do not depend on these assumptions and can therefore be evaluated more reliably (Rippe et al., 2021). In our previous study (Xue et al., 2014), the population demographic expansion of A. japonica started approximately 160,000 years ago assuming an mtCOI divergence rate of 2.4%/Myr, which was earlier than those of A. lischkeana (90,000 years ago), Atrina sp. (60,000 years ago), and A. pectinata (60,000 years ago). These results supported the following hypothesis: the East China Sea and the South China Sea have been proposed as independent marine refugia during glaciation, and Pleistocene glaciations more drastically affected species in temperate regions than in tropical regions (Hewitt, 2000;Ni et al., 2014). The expansion time of the northern group of Siphonaria japonica limpet (Gastropoda: Siphonariidae) (313,900∼682,400 years ago) was reported to be earlier than that of the southern group (104,400∼259,300 years ago) (Wang et al., 2015).

The Genetic Population Structure of the Three Atrina Species From Southern China
Since we sampled landscapes within the distributions of the three studied species along the coasts of the South China Sea, we were able to compare the population structures within the three species. We discovered significant genetic differentiation of A. pectinata between Zhanjiang and Beibu Gulf (the CT population and BH population), but this differentiation was not detected in the A. lischkeana or Atrina sp. populations in southern China. In our previous study, A. chinensis in southern China contained two lineages: the Beibu Gulf lineage and the northern South China Sea lineage (Xue et al., 2012). The Beibu Gulf is a semienclosed gulf with only one sea gate in the Qiongzhou Strait connected to the northern South China Sea. The seasonal mean current through the Qiongzhou Strait is oriented eastward from late spring through summer and westward during the rest of the year. In the Beibu Gulf, the spawning peaks of the A. pectinata species complex generally occur in the middle of May and October, and spawning troughs occur between June and July. However, we did not obtain accurate spawning periods for the three Atrina species, and we speculated that the formation of different genetic population structures among the three Atrina species may be related to their own life history characteristics (spawning period, larval diffusivity, fixed habitat, mobility, etc.). Similar population structure patterns to those of A. pectinata in southern China were also detected in Macridiscus multifarius (Bivalvia: Veneridae) (Ye et al., 2015).

CONCLUSION
The results of our study support the validity of A. japonica, A. lischkeana, and Atrina sp. as separate taxonomic units instead of as a synonyms of A. pectinata. These findings should provide useful information for developing appropriate conservation and management strategies for the A. pectinata species complex in China. Sympatric speciation illustrates how natural and sexual selection may create new species in isolation without geographic barriers. Future analyses of genome-wide diversity in these three Atrina species could provide further insight into the strength, formation and maintenance of oceanic dispersal barriers.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the NCBI GenBank under accession numbers MN817547-MN817579.

ETHICS STATEMENT
The animal study was reviewed and approved by the Animal Care Quality Assurance of the Institute of Oceanology, Chinese Academy of Sciences.

AUTHOR CONTRIBUTIONS
D-XX and TZ conceived the project. D-XX conducted the molecular experiments, data collection, analysis, and writing. H-YW assisted in data analysis and writing, provided the partial samples, assisted data interpretation, commented the results, and modified the manuscript to final version. All authors contributed to the article and approved the submitted version.