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 (nr28S rRNA: n = 44, mt16S rRNA: n = 52). Reactions were performed in 50 μl with final concentrations of 2.0 mM MgCl₂, 150 μM of each dNTP, 0.2 μM of each primer, 20 ng template DNA, 2.5 U Taq polymerase (TaKaRa; Dalian, China), and 5 μl 10 × buffer. The fragments were amplified under the following conditions: initial denaturing at 95°C for 5 min; 30 cycles of 95°C for 1 min, 58°C (nr28S rRNA), 50°C (mtCOI), 48°C (mt16S rRNA) for 1 min, and 72°C for 1 min; and a final extension at 72°C for 5 min. A negative control (no template) was included during each PCR run. The PCR products were purified using the TIANgel Midi Purification Kit (Tiangen; Beijing, China). All amplified DNA fragments were sequenced in both directions using the same PCR primers on an ABI3730 XL DNA Sequencer. Additionally, mtCOI sequences of 117 individuals collected from 6 locations along the coast of the South China Sea (Liu et al., 2011) were retrieved from GenBank (Table 1 and Figure 1), and three gene fragments for A. chinensis, A. japonica, A. exusta, A. inflata, A. rigida, and A. vexillum, and Pinna bicolor were retrieved from GenBank.

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

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 intra- and interspecific distance variations (Puillandre et al., 2012). The ABGD analysis was conducted using the online ABGD server¹ 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² 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.

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

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

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

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.

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.

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 0.000 (A. lischkeana and A. japonica) to 0.006 (Atrina sp. and A. pectinata). The pairwise genetic distances among A. lischkeana, Atrina sp., A. pectinata and A. japonica based on the mtCOI gene fragments were ranged from 0.082 (A. lischkeana and A. japonica) to 0.127 (A. lischkeana and A. pectinata). The pairwise genetic distances among A. lischkeana, Atrina sp., A. pectinata and A. japonica based on the mt16S rRNA gene fragments were ranged from 0.018 (A. lischkeana and A. japonica) to 0.059 (Atrina sp. and A. pectinata). The net genetic distances of the four species were also calculated by MEGA 6.06 as follows: A. japonica/A. lischkeana, 0.081; A. japonica/A. pectinata, 0.110; A. japonica/Atrina sp., 0.117; A. lischkeana/A. pectinata, 0.127; A. lischkeana/Atrina sp., 0.099; and A. pectinata/Atrina sp., 0.114. Applying the mtCOI divergence rate (0.7–2.4%/million years) (Liu et al., 2011), divergence events among these genealogical clades occurred approximately 3.4–18.1 million years ago during the Neogene period.

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.

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 of the three species was composed of a central haplotype that corresponded to the haplotype shared by most individuals (A. lischkeana: 61.11%, A. pectinata: 70.59%, and Atrina sp.: 70.59%) and the remaining descendant haplotypes. Within each of the three species, no distinct phylogeographic apportioning localities were observed (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.

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), and 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. lischkeana as a valid species.

Atrina sp. corresponds to the pen shells with green morphological forms (Yu et al., 2000) and with lineage 3 in Liu et al. (2011). Atrina sp. comprises rather small, thin, wedge-shaped 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.

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 (Li et al., 2015); 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).

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