A tale of two soft-shell clams: an integrative taxonomic analysis confirms Mya japonica as a valid species distinct from Mya arenaria (Bivalvia: Myidae)

A tale of two soft-shell clams: an integrative taxonomic analysis confirms Mya japonica as a... Abstract The soft-shell clam Mya arenaria Linnaeus, 1758 is a commercially important fishery resource that occurs in boreal and temperate environments in the Northern Hemisphere. Whether the soft-shell clam is a single species with a circumboreal range or a species complex also comprising Mya japonica Jay, 1857 distributed in the north Pacific has long been debated by malacologists and palaeontologists based on slight differences in shell morphology. We used an integrative taxonomic approach incorporating available Mya spp. mitochondrial COI and 16S rRNA, and nuclear 28S rRNA gene sequences, as well as spermatozoan and shell morphological characters to test the validity of M. japonica and examine the range of soft-shell clam distribution. Although differences in shell morphology were minor, the results from tree topologies, pairwise uncorrected p-distances, Automatic Barcode Gap Discovery (ABGD) and spermatozoan ultramorphological data confirm the validity of M. japonica in both its endemic region in the northwest Pacific, and as here newly reported introduced populations in British Columbia in the northeast Pacific, and show that M. arenaria is distributed in the northeast Pacific, North Atlantic, Barents Sea (Arctic Ocean) and Mediterranean. We estimate these two closely related sister species diverged 4.1–12.5 Myr during early Pliocene to late Miocene, which is consistent with current evolutionary theory regarding M. arenaria. In addition, ABGD indicated the congener Mya truncata Linnaeus, 1758 may represent a species complex, but additional evidence is still needed to clarify its taxonomic status. biogeography, cryptic species, DNA barcoding, molecular phylogeny, Mollusca, spermatozoan ultramorphology INTRODUCTION The soft-shell clam Mya arenaria Linnaeus, 1758, described from Europe, is an edible bivalve currently found in boreal and northern temperate marine and estuarine environments with a complex, but widely accepted migration history (MacNeil, 1965; Laursen, 1966; Strauch, 1972; Strasser, 1999). Mya arenaria originated in the western Pacific probably around Japan during the Miocene, spreading to the eastern Pacific later in the Miocene and then to both coasts of the Atlantic Ocean during the Pliocene (MacNeil, 1965). It was extirpated in the eastern Pacific and eastern Atlantic during the Pleistocene glaciation, surviving only in the western Atlantic with uncertain status in the Bering Sea and western Pacific (MacNeil, 1965; Carlton, 1979; Strasser, 1999). It re-invaded the eastern Atlantic during the Middle Ages and the eastern Pacific during the late 1800s, both invasions mediated by anthropogenic activities (Fujie, 1957, 1962; MacNeil, 1965; Bernard, 1979; Strasser, 1999; Cross et al., 2016). It has more recently appeared in the Mediterranean and Black Sea (Strasser, 1999; Gomoiu et al., 2002). Historical Arctic records are now generally regarded as belonging to a congener, Mya pseudoarenaria Schlesch, 1931 (Laursen, 1966). Mya arenaria is a major fishery resource in the western Atlantic (Congleton et al., 2006). However, in many parts of its introduced range M. arenaria can negatively impact native benthic communities (Streftaris & Zenetos, 2006; Crocetta & Turolla, 2011), disrupt native fisheries (Carlton, 1979) and represent a significant portion of benthic biomass (Strasser, 1999). Within the genus Mya there is much disagreement regarding evolutionary history, species distributions and taxonomy (MacNeil, 1965; Bernard, 1979; Petersen, 1999; Huber, 2010; Bouchet & Gofas, 2013; ITIS, 2017), particularly regarding the number of extant species and the relationships between extant forms and species described from fossil materials (Petersen, 1999). The most extensively studied species (both ecologically and taxonomically) is M. arenaria, but questions regarding taxonomic status and geographical range remain. Mya arenaria is reported to have numerous synonyms: M. communis Megerle von Mühlfeld, 1811; M. lata Sowerby, 1815; M. acuta Say, 1822; M. acuta mercenaria Say, 1822; M. subovata Woodward, 1833; M. subtruncata Woodward, 1833; M. alba Agassiz, 1839; M. corpulenta Conrad, 1845; M. japonica Jay, 1856; Sphenia ovoidea Carpenter, 1864; M. hemphilli Newcomb, 1874; M. declivis Pennant, 1777; M. elongata Locard, 1886; M. arenaria var. ovata Jensen, 1900; M. paternalis Matsumoto, 1930; M. arenaria oonogai Makiyama, 1935; M. japonica oonogai Makiyama, 1935; M. oonogai Makiyama, 1935; and M. arenaria corbuloides Comfort, 1938 (Coan, Scott & Bernard, 2000; Zhang, Xu & Liu, 2012). One of the most contentious questions pertains to the validity of M. japonica, which has been reported in the western Pacific (Habe, 1955, 1977; Fujie, 1957, 1962; MacNeil, 1965; Okutani, 2000), with possible relict populations or adventitious individuals occurring in the Bering Sea (MacNeil, 1965; Carlton, 1979). Soft-shell clams in the western Pacific identified as M. japonica possess rough commarginal wrinkles, relatively shorter and more rounded posterior end and more impressed pallial line in comparison to M. arenaria (Jay, 1856; MacNeil, 1965). These differences have led some malacologists to accept M. japonica as a valid species (Habe, 1955, 1977; Fujie, 1957, 1962; MacNeil, 1965; Okutani, 2000), whereas others have reported that M. arenaria and M. japonica could not be separated by the statistical analysis of morphological characteristics, and therefore, M. japonica was a junior synonym of M. arenaria (Nagao & Inoue, 1941; Bernard, 1979, 1983; Coan et al., 2000; Huber, 2010; Zhang et al., 2012). Confounding these taxonomic uncertainties are conflicting descriptions of important morphological characters, such as the chondrophore in the left valve (MacNeil, 1965; Bernard, 1979). Although several studies have examined genetic relationships between North American and European populations of M. arenaria (Strasser & Barber, 2009; Layton, Martel & Hebert, 2014; Barco et al., 2016; Cross et al., 2016; Lasota et al., 2016), none have conducted a rigorous genetic analysis of M. arenaria s.l. from across its entire range (i.e. including the western Pacific). Recent molecular analyses have identified varying levels of cryptic diversity in several cosmopolitan and wide-ranging molluscs. For example, many temperate North American populations of the tellinid clam Macoma balthica Linnaeus, 1758 are now identified as Macoma petalum Valenciennes, 1821 (Väinölä, 2003; Nikula, Strelkov & Väinölä, 2007); the widespread African freshwater oyster Etheria elliptica Lamarck, 1807 consists of at least three distinct sibling species (Elderkin et al., 2016); and the aeolidiid nudibranchs Aeolidia papillosa Linnaeus, 1761 and Spurilla neapolitana Delle Chiaje, 1841 now comprise complexes of four and five sibling species, respectively (Carmona et al., 2014; Kienberger et al., 2016). Molecular analyses can be particularly powerful when combined with morphological and other biological data (e.g. spermatozoan ultramorphology; Carstensen et al., 2009) within a broader integrative taxonomic framework (Padial et al., 2010; Schlick-Steiner et al., 2010). Spermatozoan ultramorphology appears to be highly conserved and has been used to resolve taxonomic and phylogenetic questions at the species level in many bivalve genera, including Solen (Hodgson, Devilliers & Bernard, 1987), Donax (Hodgson, Bernard & Vanderhorst, 1990; Introini, Passos & Recco-Pimentel, 2013), Crassostrea (Gwo, Liou & Cheng, 1996), Bathymodiolus (Eckelbarger & Young, 1999), Nutricola (Geraghty, Russell & Dollahon, 2008) and Brachidontes (Torroglosa & Gimenez, 2015). Fine organization of soft-shell clam spermatozoa from the western Pacific was studied in M. japonica (Drozdov, Sharina & Tyurin, 2009) and M. arenaria oonogai (Kim, Chung & Park, 2011). In both studies, spermatozoa were classified as ect-aquasperm with a conical curved head consisting of an acrosomal complex and nucleus, a midpiece with four mitochondria and a flagellum with a 9 × 2 + 2 organization of microtubules. However, no previous studies have described the morphology of M. arenaria spermatozoa from the greater Atlantic or preformed range-wide comparisons. In this study, we report the first rigorous range-wide investigation of the soft-shell clam, using an integrative taxonomic approach based on mitochondrial cytochrome c oxidase subunit I (COI), mitochondrial 16S rRNA (16S) and nuclear 28S rRNA (28S) gene data, and quantitative sperm ultramorphology and qualitative shell morphological comparisons. We also report the first detailed description of spermatozoan ultramorphology of M. arenaria s.s. and provide updated distributional data. MATERIAL AND METHODS Specimen collection For DNA analyses, we collected soft-shell clams from Yellow Sea, Qingdao, China (n = 8); Barents Sea, Murmansk, Russia (n = 3); Vostok Bay, Sea of Japan, Russia (n = 1); and Chesapeake Bay, Maryland, USA (n = 10) (Fig. 1). We also collected a congener, Mya uzenensis Nomura et Zinbo, 1937 from Vostok Bay (n = 1). The authors follow Lutaenko & Noseworthy (2012) and accept M. uzenensis as a valid species and not a junior synonym of M. pseudoarenaria. Voucher specimens were deposited at the Marine Biological Museum, Chinese Academy of Sciences (MBM), Qingdao, China, and the Smithsonian Institution National Museum of Natural History (USNM), Washington, DC, USA (Table 1). For spermatozoan ultramorphological analyses, three males in each of the Yellow Sea, Barents Sea and Sea of Japan sampling sites were collected at pre-spawning and spawning stages of the reproductive cycle (for staging details, see Porter, 1974; Dzyuba & Maslennikova, 1987). Voucher specimens were deposited at the Zoological Museum, Far Eastern Federal University, Vladivostok, Russia (ZMFU). In addition to the samples mentioned above, we examined 135 lots of soft-shell clams (photographing a sub-set of samples), which included individuals previously identified as M. arenaria and M. japonica, from historical MBM and ZMFU collections to identify possible morphological differences between soft-shell clams collected in the western Pacific and northwestern Asia/northern Europe (e.g. Barents Sea, White Sea, North Sea, Baltic Sea; Fig. 2). Table 1. Sampling information of specimens collected for this study and GenBank accession numbers of their sequences Species Voucher Location GenBank numbers COI 16S 28S Mya arenaria MBM229022 Murmansk, Barents Sea, Russia (Arctic) KX534201 KX534176 KX534188 Mya arenaria MBM229023 Murmansk, Barents Sea, Russia (Arctic) KX534202 KX534177 KX534189 Mya arenaria MBM229024 Murmansk, Barents Sea, Russia (Arctic) KX534203 KX534178 KX534190 Mya arenaria USNM:IZ:1286707 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959460 KT959520 – Mya arenaria USNM:IZ:1286708 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959399 KT959490 – Mya arenaria USNM:IZ:1286709 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959418 KT959498 – Mya arenaria USNM:IZ:1286710 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959394 KT959487 – Mya arenaria USNM:IZ:1286892 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906055 – – Mya arenaria USNM:IZ:1286895 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906064 – – Mya arenaria USNM:IZ:1286896 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905943 – – Mya arenaria USNM:IZ:1286897 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905735 – – Mya arenaria USNM:IZ:1287437 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906072 – – Mya arenaria USNM:IZ:1287540 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906038 – – Mya japonica MBM229013 Qingdao, Yellow Sea, China (NW Pacific) KX534192 KX534167 KX534180 Mya japonica MBM229014 Qingdao, Yellow Sea, China (NW Pacific) KX534193 KX534168 – Mya japonica MBM229015 Qingdao, Yellow Sea, China (NW Pacific) KX534194 KX534169 KX534181 Mya japonica MBM229016 Qingdao, Yellow Sea, China (NW Pacific) KX534195 KX534170 KX534182 Mya japonica MBM229017 Qingdao, Yellow Sea, China (NW Pacific) KX534196 KX534171 KX534183 Mya japonica MBM229018 Qingdao, Yellow Sea, China (NW Pacific) KX534197 KX534172 KX534184 Mya japonica MBM229019 Qingdao, Yellow Sea, China (NW Pacific) KX534198 KX534173 KX534185 Mya japonica MBM229020 Qingdao, Yellow Sea, China (NW Pacific) KX534199 KX534174 KX534186 Mya japonica MBM229021 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534200 KX534175 KX534187 Mya uzenensis MBM229025 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534204 KX534179 KX534191 Species Voucher Location GenBank numbers COI 16S 28S Mya arenaria MBM229022 Murmansk, Barents Sea, Russia (Arctic) KX534201 KX534176 KX534188 Mya arenaria MBM229023 Murmansk, Barents Sea, Russia (Arctic) KX534202 KX534177 KX534189 Mya arenaria MBM229024 Murmansk, Barents Sea, Russia (Arctic) KX534203 KX534178 KX534190 Mya arenaria USNM:IZ:1286707 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959460 KT959520 – Mya arenaria USNM:IZ:1286708 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959399 KT959490 – Mya arenaria USNM:IZ:1286709 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959418 KT959498 – Mya arenaria USNM:IZ:1286710 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959394 KT959487 – Mya arenaria USNM:IZ:1286892 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906055 – – Mya arenaria USNM:IZ:1286895 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906064 – – Mya arenaria USNM:IZ:1286896 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905943 – – Mya arenaria USNM:IZ:1286897 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905735 – – Mya arenaria USNM:IZ:1287437 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906072 – – Mya arenaria USNM:IZ:1287540 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906038 – – Mya japonica MBM229013 Qingdao, Yellow Sea, China (NW Pacific) KX534192 KX534167 KX534180 Mya japonica MBM229014 Qingdao, Yellow Sea, China (NW Pacific) KX534193 KX534168 – Mya japonica MBM229015 Qingdao, Yellow Sea, China (NW Pacific) KX534194 KX534169 KX534181 Mya japonica MBM229016 Qingdao, Yellow Sea, China (NW Pacific) KX534195 KX534170 KX534182 Mya japonica MBM229017 Qingdao, Yellow Sea, China (NW Pacific) KX534196 KX534171 KX534183 Mya japonica MBM229018 Qingdao, Yellow Sea, China (NW Pacific) KX534197 KX534172 KX534184 Mya japonica MBM229019 Qingdao, Yellow Sea, China (NW Pacific) KX534198 KX534173 KX534185 Mya japonica MBM229020 Qingdao, Yellow Sea, China (NW Pacific) KX534199 KX534174 KX534186 Mya japonica MBM229021 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534200 KX534175 KX534187 Mya uzenensis MBM229025 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534204 KX534179 KX534191 View Large Table 1. Sampling information of specimens collected for this study and GenBank accession numbers of their sequences Species Voucher Location GenBank numbers COI 16S 28S Mya arenaria MBM229022 Murmansk, Barents Sea, Russia (Arctic) KX534201 KX534176 KX534188 Mya arenaria MBM229023 Murmansk, Barents Sea, Russia (Arctic) KX534202 KX534177 KX534189 Mya arenaria MBM229024 Murmansk, Barents Sea, Russia (Arctic) KX534203 KX534178 KX534190 Mya arenaria USNM:IZ:1286707 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959460 KT959520 – Mya arenaria USNM:IZ:1286708 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959399 KT959490 – Mya arenaria USNM:IZ:1286709 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959418 KT959498 – Mya arenaria USNM:IZ:1286710 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959394 KT959487 – Mya arenaria USNM:IZ:1286892 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906055 – – Mya arenaria USNM:IZ:1286895 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906064 – – Mya arenaria USNM:IZ:1286896 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905943 – – Mya arenaria USNM:IZ:1286897 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905735 – – Mya arenaria USNM:IZ:1287437 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906072 – – Mya arenaria USNM:IZ:1287540 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906038 – – Mya japonica MBM229013 Qingdao, Yellow Sea, China (NW Pacific) KX534192 KX534167 KX534180 Mya japonica MBM229014 Qingdao, Yellow Sea, China (NW Pacific) KX534193 KX534168 – Mya japonica MBM229015 Qingdao, Yellow Sea, China (NW Pacific) KX534194 KX534169 KX534181 Mya japonica MBM229016 Qingdao, Yellow Sea, China (NW Pacific) KX534195 KX534170 KX534182 Mya japonica MBM229017 Qingdao, Yellow Sea, China (NW Pacific) KX534196 KX534171 KX534183 Mya japonica MBM229018 Qingdao, Yellow Sea, China (NW Pacific) KX534197 KX534172 KX534184 Mya japonica MBM229019 Qingdao, Yellow Sea, China (NW Pacific) KX534198 KX534173 KX534185 Mya japonica MBM229020 Qingdao, Yellow Sea, China (NW Pacific) KX534199 KX534174 KX534186 Mya japonica MBM229021 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534200 KX534175 KX534187 Mya uzenensis MBM229025 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534204 KX534179 KX534191 Species Voucher Location GenBank numbers COI 16S 28S Mya arenaria MBM229022 Murmansk, Barents Sea, Russia (Arctic) KX534201 KX534176 KX534188 Mya arenaria MBM229023 Murmansk, Barents Sea, Russia (Arctic) KX534202 KX534177 KX534189 Mya arenaria MBM229024 Murmansk, Barents Sea, Russia (Arctic) KX534203 KX534178 KX534190 Mya arenaria USNM:IZ:1286707 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959460 KT959520 – Mya arenaria USNM:IZ:1286708 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959399 KT959490 – Mya arenaria USNM:IZ:1286709 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959418 KT959498 – Mya arenaria USNM:IZ:1286710 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959394 KT959487 – Mya arenaria USNM:IZ:1286892 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906055 – – Mya arenaria USNM:IZ:1286895 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906064 – – Mya arenaria USNM:IZ:1286896 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905943 – – Mya arenaria USNM:IZ:1286897 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905735 – – Mya arenaria USNM:IZ:1287437 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906072 – – Mya arenaria USNM:IZ:1287540 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906038 – – Mya japonica MBM229013 Qingdao, Yellow Sea, China (NW Pacific) KX534192 KX534167 KX534180 Mya japonica MBM229014 Qingdao, Yellow Sea, China (NW Pacific) KX534193 KX534168 – Mya japonica MBM229015 Qingdao, Yellow Sea, China (NW Pacific) KX534194 KX534169 KX534181 Mya japonica MBM229016 Qingdao, Yellow Sea, China (NW Pacific) KX534195 KX534170 KX534182 Mya japonica MBM229017 Qingdao, Yellow Sea, China (NW Pacific) KX534196 KX534171 KX534183 Mya japonica MBM229018 Qingdao, Yellow Sea, China (NW Pacific) KX534197 KX534172 KX534184 Mya japonica MBM229019 Qingdao, Yellow Sea, China (NW Pacific) KX534198 KX534173 KX534185 Mya japonica MBM229020 Qingdao, Yellow Sea, China (NW Pacific) KX534199 KX534174 KX534186 Mya japonica MBM229021 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534200 KX534175 KX534187 Mya uzenensis MBM229025 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534204 KX534179 KX534191 View Large Figure 1. View largeDownload slide Locations of Mya japonica (red triangles), Mya arenaria (red dots) and Mya uzenensis (red diamond) collected in this study, and other known localities of Mya japonica (yellow triangles), Mya arenaria (yellow dots) and Mya truncata (yellow squares) of which sequences were used in genetic analyses. Figure 1. View largeDownload slide Locations of Mya japonica (red triangles), Mya arenaria (red dots) and Mya uzenensis (red diamond) collected in this study, and other known localities of Mya japonica (yellow triangles), Mya arenaria (yellow dots) and Mya truncata (yellow squares) of which sequences were used in genetic analyses. Figure 2. View largeDownload slide Mya shell morphology. A–E, Mya japonica Jay, 1857, ZMFU no. 10069/Bv-514, sampled from the open part of Peter the Great Bay, Bolshoy Pelis Island, Russia. F–J, Mya japonica Jay, 1857, ZMFU no. 23301/Bv-3785, sampled from semi-enclosed part of Peter the Great Bay, Vityaz Bay, Russia. K–O, Mya japonica Jay, 1857, MBM229023, sampled from Qingdao, Yellow Sea, China. P–T, Mya arenaria Linnaeus, 1758, MBM229013, sampled from Murmansk, Barents Sea, Russia. Figure 2. View largeDownload slide Mya shell morphology. A–E, Mya japonica Jay, 1857, ZMFU no. 10069/Bv-514, sampled from the open part of Peter the Great Bay, Bolshoy Pelis Island, Russia. F–J, Mya japonica Jay, 1857, ZMFU no. 23301/Bv-3785, sampled from semi-enclosed part of Peter the Great Bay, Vityaz Bay, Russia. K–O, Mya japonica Jay, 1857, MBM229023, sampled from Qingdao, Yellow Sea, China. P–T, Mya arenaria Linnaeus, 1758, MBM229013, sampled from Murmansk, Barents Sea, Russia. Prior to analyses, all specimens (field collected and museum vouchers) were identified based on shell morphological differences of these two putative species (M. arenaria or M. japonica) recognized by early malacologists (Habe, 1955, 1977; Fujie, 1957, 1962; MacNeil, 1965). Specifically, we examined shell shape, noting profile of the posterior end, shell texture (regularity and depth of commarginal growth lines), the pallial line, and the shape, size, and orientation of the chondrophore in the left valve. Soft-shell clams possessing an elongated shell shape with even commarginal growth lines and lightly impressed pallial line were classified as M. arenaria, whereas soft-shell clams with a more truncate shell shape with rough commarginal wrinkles and more deeply impressed pallial line were classified as M. japonica. Molecular work DNA extraction, amplification and sequencing Thirteen specimens of M. arenaria, nine specimens of M. japonica and one specimen of M. uzenensis were collected for molecular analyses. We extracted total DNA from the samples using the Marine Animal Genomic DNA Extraction Kit (Tiangen Biotech, Beijing, China) according to manufacturer’s instructions. North American samples were sequenced following Aguilar et al. (2016) using the Primer sets jgLCO1490/jgHCO2198 (COI) and 16Sar-L/16Sb (16S). For Eurasian samples, fragments of other specimens of the mitochondrial COI gene were amplified using primers of Lco1490: 5′-GGT CAA CAA ATC ATA AAG ATA TTG G-3′, Hco2198: 5′-TAA ACT TCA GGG TGA CCA AAA AAT CA-3′ (Folmer et al., 1994). Fragments of the mitochondrial 16S rRNA gene were amplified using 16Sa: 5′-CGC CTG TTT ATC AAA AAC AT-3′ (Xiong & Kocher, 1991), 16Sb: 5′-CTC CGG TTT GAA CTC AGA TCA-3′ (Edgecombe, Giribet & Wheeler, 2002). The nuclear 28S rRNA gene was amplified using 28Sa: 5′-GAC CCG TCT TGA AAC ACG GA-3′, 28Sb: 5′-TCG GAA GGA ACC AGC TAC-3′ (Whiting et al., 1997). The 25-µL reaction for the amplification contained 12 µL 2×Es TaqMasterMix (CWBio Co., Ltd, Beijing, China), 2 µL of template DNA (50 ng/µL), 1 µL of each primer (10 M) and 9 µL dH2O. Amplification reactions were conducted in a T100 Thermal Cycler (Bio-Rad Laboratories, Inc.) with the following thermal profile: 95 °C for 3 min; 35 cycles of 95 °C for 30 s, 50–55 °C for 45 s, 72 °C for 1 min, and a further 10-min elongation at 72 °C. The PCR products were purified and sequenced by BGI Tech Solutions Co., Ltd. Species delimitation analyses To discriminate species, we used an integrative taxonomic approach including pairwise uncorrected p-distances, phylogenetic trees, Automatic Barcode Gap Discovery (ABGD) and spermatozoan ultramorphological and qualitative shell morphological data. In addition to sequences generated from our field collections, we included all Mya COI (> 500 bp), 16S and 28S gene sequences available in GenBank and BOLD (Ratnasingham & Hebert, 2007) in the generation of p-distances and in the ABGD analyses, most of which were associated with published articles (Supporting Information, Table S1). In total, 175 COI sequences (including four Potamocorbula sequences, which were originally submitted under the name Corbula amurensis in GenBank as an outgroup) were used for ABGD and p-distance analyses, and 26 16S sequences (including four outgroup Potamocorbula sequences) were also used for p-distance analyses. To improve clarity, sequences without associated publications or location information were omitted (n = 49) from phylogenetic trees. The ABGD (Puillandre et al., 2012) settings were as follows: steps = 10, X = 1, Pmin = 0.001, Pmax = 0.1, Nb bins = 20, and Jukes-Cantor (JC69) and Kimura (K80) TS/TV. The p-distances within and among each species grouping were estimated with MEGA v.6 (Tamura et al., 2013). Phylogenetic analyses Sequences of each region were aligned independently using MAFFT v.7 (Katoh & Standley, 2013) with the G-INS-i and Q-INS-i algorithms for the protein-coding and ribosomal regions. The alignments were further checked by eye. A combined analysis was conducted with concatenated sequences from the three genes. Three genes from the same individual were concatenated together as a sequence in SequenceMatrix v 1.8 (Vaidya, Lohman & Meier, 2011). Sequences from Potamocorbula amurensis (Schrenck, 1861) in the same superfamily, Myoidea, were selected as an outgroup. The best-fitted evolutionary models were selected by Akaike information criterion as implemented in jModeltest2 (Darriba et al., 2012). Maximum likelihood (ML) analysis was carried out using PhyML-3.1 (Guindon et al., 2010). Node support came from a majority-rule consensus tree of 1000 bootstrap replicates. For the ML bootstraps, we consider values < 70 as low, 70‒95 as moderate and ≥ 95 as high, following Hillis & Bull (1993). Bayesian inference (BI) analysis was carried out using MrBayes v3.2.6 (Ronquist et al., 2012) in CIPRES Science Gateway. We estimated posterior probability (PP) using four chains running 10000000 generations sampling every 100 generations. The first 25% of sampled trees were considered burn-in trees, and Tracer v1.6 was used to confirm that post-burn-in trees yielded an effective sample size of > 200 for all parameters. All remaining trees were used to calculate PPs using a 50% majority-rule consensus. For the Bayesian PPs, we considered values < 94% as low and ≥ 95% as high following Alfaro, Zoller & Lutzoni (2003). Divergence time estimates We used a Bayesian Markov chain Monte Carlo method implemented in Beast version 2.4.0 (Bouckaert et al., 2014) to estimate the evolutionary ages of internal nodes in the tree topology derived from the combined phylogenetic analyses. The time to most recent common ancestor (T), that is the time of divergence, was estimated as T = K/2r, where K is the K2P genetic distance, r is the nucleotide substitution rate. Substitution rates (r) were estimated as percentage per lineage per million years (Myr), so they equal one half the pairwise sequence divergence per Myr, that is divergence rate. Our estimates of T are based on the mitochondrial COI divergence between cognate species of Arcidae isolated across the Isthmus of Panama at a rate of 0.7–1.2%/Myr (Marko, 2002); thus, the substitution rate (r) equals 0.35–0.6%/Myr. A relaxed, uncorrelated lognormal clock with a general time reversible (GTR) substitution model was inferred for each partition, and a Yule speciation process was assumed for the tree prior. The analysis was run for 100 million generations, sampled every 1000 generations and a burn-in of 25%. Convergence and mixing were assessed using Tracer v. 1.6 by examining log-likelihood values across generations and ensuring that all remaining trees yielded an effective sample size for all parameters. Results were visualized using Tracer v. 1.6 and FigTree v. 1.4.2. Spermatozoan ultramorphology Transmission electron microscopes preparation Small pieces of gonadal tissues were pre-fixed in 2.5% (v/v) glutaraldehyde in 0.1 M cacodylate buffer (pH 7.5), with the addition of NaCl to provide appropriate tonicity of sea water, for 2 h at 4 °С. The specimens were then rinsed several times in the same buffer and held at 4 °C until analysis at the National Scientific Center of Marine Biology (Vladivostok, Russia). In the laboratory, gonadal tissue samples were post-fixed in 2% (w/v) osmium tetroxide in 0.1 M cacodylate buffer (pH 7.5) for 2 h in the dark at room temperature. The specimens were washed in distilled water and then dehydrated using a graded ethanol series and acetone. The dehydrated specimens were embedded in Araldite-EMbed-812 resin. Ultrathin sections (c. 50 nm) were obtained using a Leica UC6 ultramicrotome equipped with a diamond knife. The sections were stained with 2% aqueous uranyl acetate and Reynolds’ lead citrate (Reynolds, 1963) and were observed under Zeiss Libra 120 and Zeiss Libra 200FE transmission electron microscopes (Far East Center of Electron Microscopy, National Scientific Center of Marine Biology FEB RAS). Measurements of spermatozoa To investigate possible differences in spermatozoan ultramorphology among soft-shell clams, we examined three males (five spermatozoa per male) from each species/geographical groupings: M. arenaria from the Barents Sea, M. japonica from the Sea of Japan and M. japonica from the Yellow Sea. Only those spermatozoa in which the following compartments could be simultaneously observed in a longitudinal projection were measured: (1) acrosomal complex with cylindrical subacrosomal space, from the base of anterior nuclear fossa to the apex of acrosomal vesicle; (2) cone-shaped nucleus with the gradual dilatation from the apical to the basal part without any narrowing; and (3) midpiece with two centrioles. Because of the curved shape of the nucleus, we measured the maximal Feret’s diameter to estimate the nuclear length. We measured the length of the acrosomal complex from the invagination of the anterior nuclear fossa to the terminal point at the apical part of the acrosomal vesicle, and the width of the acrosomal complex at its base. All measurements were calculated with ImageJ version 1.49b (Schneider, Rasband & Eliceiri, 2012). To compare the size and shape of the acrosomal complex, we calculated volumes and length to width ratios. Based on the close-to-conical profile of the acrosomal complexes in the sections, we used the following formula to estimate volume: V = 1/3πR2H, where R is ½ of the width of the acrosomal complex and H is the length of the acrosomal complex. Statistical analyses We used one-way analysis of variance (ANOVA) to test for significant differences in mean nucleus length, mean acrosomal volume and mean acrosomal length to width ratio among the three species/geographical groupings: M. arenaria from the Barents Sea, M. japonica from the Sea of Japan and M. japonica from the Yellow Sea. Effect size (η2) for ANOVA was calculated using the corresponding function of R add-on package ‘lsr’ version 0.5 (Navarro, 2014). Multiple comparisons among geographical areas were made using Tukey’s honestly significant difference (Tukey’s HSD) tests. Prior to analysis, all data sets were tested for normality with Shapiro-Wilk tests and for homogeneity of variance using Levene’s tests and were determined to meet these assumptions of parametric analyses. Differences were considered statistically significant at P < 0.05 for all tests. We performed all statistical analyses using R version 3.2.4 (R Core Team, 2015) running on a Debian GNU/Linux workstation. RESULTS Morphology Shell characteristics and morphological identifications Based on our visual examinations, the samples from the western Pacific morphologically corresponded well with M. japonica and the northeast Pacific/Atlantic/Arctic soft-shell calms were morphologically identical with M. arenaria (Fig. 2). Detailed statistical analyses of morphological characteristics of soft-shell clams classified as M. arenaria and M. japonica have been reported previously (see Bernard, 1979), but here we describe minor differences in shell characteristics between them. In general, M. arenaria possessed a more elongated shell, with a longer and thinner posterior end and regular thin commarginal growth lines (Fig. 2P–T). Mya japonica was higher, with rough commarginal wrinkles and a relatively shorter and more rounded posterior end (Fig. 2A–O). In addition, the pallial line of M. arenaria was less impressed than in M. japonica. However, there was considerable variation in the shape, relative size and orientation of the left valve chondrophore among all the clams examined and not useful for discriminating between species. Two distinct morphological forms of M. japonica were observed in Peter the Great Bay. A thin, ‘typical’ form of M. japonica, which was elongated with regular growth lines, occurring mostly in semi-enclosed bays and other protected seashore areas, often with decreased salinity and a low energy environment (Fig. 2F–J). Conversely, a ‘thick-shelled’ form with irregular growth lines occurred in open, wave-influenced bays (Fig. 2A–E). Intermediate specimens were also found. Comparison of spermatozoan ultramorphology All spermatozoa examined in the present study were ect-aquasperm, a presumed primitive sperm type often associated with aquatic invertebrates with external fertilizing, free-swimming gametes (Eckelbarger & Young, 1999). Spermatozoa of all the specimens examined consisted of an acrosomal complex, nucleus, midpiece and flagellum. The data on dimensions of nuclei and acrosomal complexes are summarized in Table 2. The sperm were all rotationally asymmetrical due to the curvature of the nucleus or the position of the acrosomal complex (Fig. 3A–C). The midpieces of all spermatozoa examined were similar, containing four spherical mitochondria (Fig. 3G–I), and orthogonally arranged proximal and distal centrioles (Fig. 3J–L). In addition, the flagellum of all spermatozoa possessed a typically organized axoneme (9 × 2 + 2), posterior to the distal centriole. Table 2. Descriptive statistics of spermatozoan compartment dimensions in Mya arenaria sampled from the Barents Sea and in Mya japonica sampled from the Sea of Japan and the Yellow Sea (mean ± SD) Mya arenaria from the Barents Sea Mya japonica from the Sea of Japan Mya japonica from the Yellow Sea Nucleus length (µm) 1.98 ± 0.16 2.97 ± 0.13 3.08 ± 0.18 Acrosomal complex length (µm) 0.74 ± 0.06 0.80 ± 0.05 0.82 ± 0.04 Acrosomal complex width (µm) 0.61 ± 0.05 0.55 ± 0.04 0.57 ± 0.03 Acrosomal complex length to width ratio 1.22 ± 0.11 1.46 ± 0.08 1.46 ± 0.09 Acrosomal complex volume (µm3) 0.073 ± 0.015 0.064 ± 0.012 0.069 ± 0.010 Mya arenaria from the Barents Sea Mya japonica from the Sea of Japan Mya japonica from the Yellow Sea Nucleus length (µm) 1.98 ± 0.16 2.97 ± 0.13 3.08 ± 0.18 Acrosomal complex length (µm) 0.74 ± 0.06 0.80 ± 0.05 0.82 ± 0.04 Acrosomal complex width (µm) 0.61 ± 0.05 0.55 ± 0.04 0.57 ± 0.03 Acrosomal complex length to width ratio 1.22 ± 0.11 1.46 ± 0.08 1.46 ± 0.09 Acrosomal complex volume (µm3) 0.073 ± 0.015 0.064 ± 0.012 0.069 ± 0.010 Number of samples: 15 spermatozoa in each sampling site. View Large Table 2. Descriptive statistics of spermatozoan compartment dimensions in Mya arenaria sampled from the Barents Sea and in Mya japonica sampled from the Sea of Japan and the Yellow Sea (mean ± SD) Mya arenaria from the Barents Sea Mya japonica from the Sea of Japan Mya japonica from the Yellow Sea Nucleus length (µm) 1.98 ± 0.16 2.97 ± 0.13 3.08 ± 0.18 Acrosomal complex length (µm) 0.74 ± 0.06 0.80 ± 0.05 0.82 ± 0.04 Acrosomal complex width (µm) 0.61 ± 0.05 0.55 ± 0.04 0.57 ± 0.03 Acrosomal complex length to width ratio 1.22 ± 0.11 1.46 ± 0.08 1.46 ± 0.09 Acrosomal complex volume (µm3) 0.073 ± 0.015 0.064 ± 0.012 0.069 ± 0.010 Mya arenaria from the Barents Sea Mya japonica from the Sea of Japan Mya japonica from the Yellow Sea Nucleus length (µm) 1.98 ± 0.16 2.97 ± 0.13 3.08 ± 0.18 Acrosomal complex length (µm) 0.74 ± 0.06 0.80 ± 0.05 0.82 ± 0.04 Acrosomal complex width (µm) 0.61 ± 0.05 0.55 ± 0.04 0.57 ± 0.03 Acrosomal complex length to width ratio 1.22 ± 0.11 1.46 ± 0.08 1.46 ± 0.09 Acrosomal complex volume (µm3) 0.073 ± 0.015 0.064 ± 0.012 0.069 ± 0.010 Number of samples: 15 spermatozoa in each sampling site. View Large Figure 3. View largeDownload slide Sperm fine morphology in the studied Mya species/populations. A, D, G, J, Mya arenaria sampled from the Barents Sea. B, E, H, K, Mya japonica sampled from the Sea of Japan. C, F, I, L, Mya japonica sampled from the Yellow Sea. A–C, longitudinal section of spermatozoon consisted of acrosomal complex (ac), nucleus (n) with anterior (anf) and posterior nuclear fossae (pnf), midpiece (mp) and flagellum (f), scale bar = 1 µm. D–F, acrosomal complex included acrosomal vesicle (av) with an inner moderate electron-dense layer (asterisk) and outer electron-dense layer (arrow) and subacrosomal space (sas) filled with subacrosomal material, scale bar = 0.3 µm. G–I, transversal section of midpiece showing four mitochondria (m) around distal centriole (dc), scale bar = 0.5 µm. J–L, longitudinal section of midpiece. Proximal (pc) and distal centrioles (dc) located orthogonally to each other, scale bar = 0.5 µm. Figure 3. View largeDownload slide Sperm fine morphology in the studied Mya species/populations. A, D, G, J, Mya arenaria sampled from the Barents Sea. B, E, H, K, Mya japonica sampled from the Sea of Japan. C, F, I, L, Mya japonica sampled from the Yellow Sea. A–C, longitudinal section of spermatozoon consisted of acrosomal complex (ac), nucleus (n) with anterior (anf) and posterior nuclear fossae (pnf), midpiece (mp) and flagellum (f), scale bar = 1 µm. D–F, acrosomal complex included acrosomal vesicle (av) with an inner moderate electron-dense layer (asterisk) and outer electron-dense layer (arrow) and subacrosomal space (sas) filled with subacrosomal material, scale bar = 0.3 µm. G–I, transversal section of midpiece showing four mitochondria (m) around distal centriole (dc), scale bar = 0.5 µm. J–L, longitudinal section of midpiece. Proximal (pc) and distal centrioles (dc) located orthogonally to each other, scale bar = 0.5 µm. Nuclei were curved in all three geographical groups, but in M. arenaria from the Barents Sea, the curvature of the nuclei was less obvious than that in M. japonica from the Sea of Japan and the Yellow Sea (Fig. 3A–C). All nuclei possessed electron-dense chromatin with irregularly shaped electron-lucent lacunae. Anterior portions of the nuclei were slightly invaginated and formed an anterior nuclear fossa that was located subterminally (in M. arenaria) or apically (in M. japonica; Fig. 3A–C). Posterior portions had a well-defined invagination (posterior nuclear fossa) where a proximal centriole was located, and indentations were in the regions of contact of the nucleus and mitochondria (Fig. 3A–C, J–L). There were statistically significant differences in mean nuclei length (ANOVA: F(2,42) = 219.51, P < 0.05, η2 = 0.91). The Tukey’s HSD test indicated that M. arenaria (Barents Sea) had significantly shorter nuclei compared to M. japonica from the Sea of Japan and Yellow Sea (P < 0.05; Fig. 4A). No statistically significant difference in nuclei length between M. japonica collected from the Sea of Japan and Yellow Sea was found (P = 0.165). Figure 4. View largeDownload slide Whisker plots showing results of Tukey’s honestly significant difference test for difference in mean nucleus lengths (A) and acrosomal L/W ratios (B) among groups of Mya. Black dots are the difference in pairwise comparisons between sampling locations of the means and the whiskers are 95% confidence intervals of these differences. Figure 4. View largeDownload slide Whisker plots showing results of Tukey’s honestly significant difference test for difference in mean nucleus lengths (A) and acrosomal L/W ratios (B) among groups of Mya. Black dots are the difference in pairwise comparisons between sampling locations of the means and the whiskers are 95% confidence intervals of these differences. In all spermatozoa examined, the acrosomal complexes were cone-shaped and consisted of a membrane-bounded acrosomal vesicle and subacrosomal space (Fig. 3D–F). Within the acrosomal vesicle, two regions with different electron density were distinguished: the inner with moderate electron density and the outer with high electron density (Fig. 3D–F). The acrosomal vesicle had a deep posterior invagination that extended to its apical portion and formed a subacrosomal space filled with flake-like subacrosomal material of low electron density. The distance between the acrosomal complex and the nucleus was narrow in M. arenaria (Barents Sea), but wider M. japonica (Sea of Japan and Yellow Sea; Fig. 3D–F). This is most likely the reason why the acrosomal complex lengths were generally larger relative to width in M. japonica (Sea of Japan and Yellow Sea) than M. arenaria (Barents Sea; Table 2). There was a statistically significant difference in mean acrosomal length to width ratios (ANOVA: F(2,42) = 30.83, P < 0.05, η2 = 0.60). A Tukey’s HSD test indicated that acrosomal complexes of the spermatozoa in M. arenaria (Barents Sea) were shorter relative to width than those in M. japonica (Sea of Japan and Yellow Sea; Table 2; P < 0.05), while no statistically significant difference in length to width ratios between Mya from the Sea of Japan and Yellow Sea was found (P = 1; Fig. 4B). There were no significant differences in mean acrosomal complex volumes among the different geographical regions (ANOVA: F(2,42) = 1.94, P = 0.156, η2 = 0.08). Molecular analyses Barcoding analyses In total, 52 new sequences were obtained and deposited in GenBank (Table 1). The alignment data set comprised 637 nucleotide positions for COI, 496 positions for 16S and 297 positions for 28S. The concatenated data set yielded a sequence alignment of 1430 positions. The ABGD analysis based on the COI alignment produced five distinct groupings, with P values ranging from 0.022 to 0.100 based on both the Jukes-Cantor (JC69) and Kimura (K80) TS/TV models (Table 3). These groupings corresponded with the five clades/subclades in the COI trees (Table 3; Fig. 5B–D). The COI p-distances within the former four groups were low (0–2.39%), while the p-distances among the five groups ranged from 7.17 to 19.12% (Table 3; Supporting Information, Table S2). The minimum interspecific p-distance among these groups was much higher than the maximum one within ABGD Group 1 (7.17 vs. 2.39%; Table 3), which indicated all the western Pacific soft-shell clams (those morphologically identified M. japonica by the authors and those classified as M. arenaria – the valid species name at the time of identification), plus two individuals from British Columbia (DFO097-11 + RBCMI219-14) represent a single species (i.e. M. japonica). After reassigning all western Pacific M. arenaria and DFO097-11 + RBCMI219-14 as M. japonica, the COI intraspecific p-distances within Mya were 0–7.57%, while the interspecific p-distances among Mya species ranged from 11.16 to 19.12% (Fig. 5A; Table 3; Supporting Information, Table S2), indicating a barcoding gap among Mya species (Fig. 5A). However, if the two Mya truncata clades (ABGD Groups 3 and 4) are considered distinct species, intraspecific p-distances decrease to 0–2.39%, while interspecific distances range from 7.17 to 19.12%. For 16S, the intraspecific p-distances ranged from 0 to 0.97% within M. arenaria and M. japonica, while the interspecific p-distances ranged from 5.3 to 13.80% (Supporting Information, Table S3), which also indicates a barcoding gap between M. arenaria and M. japonica. Table 3. Uncorrected pairwise distances at COI within and between the ABGD groups Group 1 Group 2 Group 3 Group 4 Group 5 Group 1 All Mya japonica + Mya arenaria QWEAS092- 15, QWEAS093-15, DFO097-11, RBCMI219- 14, KP976290, KP976289, KP976288, KP976287, KP976286, KJ125420, KJ125421, JQ267795 0‒2.39% Group 2 Remaining Mya arenaria (142 sequences) 11.55–13.94% 0‒1.99% Group 3 Mya truncata KF644116, KF644129 17.53–18.33% 15.54–17.53% 0.40% Group 4 Mya truncata KF644154, KF643769, KF643675, KF643403 17.53–18.73% 17.53–19.12% 7.17–11.55% 0 Group 5 Mya uzenensis KX534204 17.53–18.33% 16.33–17.93% 11.16–11.55% 13.15% ‒ Group 1 Group 2 Group 3 Group 4 Group 5 Group 1 All Mya japonica + Mya arenaria QWEAS092- 15, QWEAS093-15, DFO097-11, RBCMI219- 14, KP976290, KP976289, KP976288, KP976287, KP976286, KJ125420, KJ125421, JQ267795 0‒2.39% Group 2 Remaining Mya arenaria (142 sequences) 11.55–13.94% 0‒1.99% Group 3 Mya truncata KF644116, KF644129 17.53–18.33% 15.54–17.53% 0.40% Group 4 Mya truncata KF644154, KF643769, KF643675, KF643403 17.53–18.73% 17.53–19.12% 7.17–11.55% 0 Group 5 Mya uzenensis KX534204 17.53–18.33% 16.33–17.93% 11.16–11.55% 13.15% ‒ View Large Table 3. Uncorrected pairwise distances at COI within and between the ABGD groups Group 1 Group 2 Group 3 Group 4 Group 5 Group 1 All Mya japonica + Mya arenaria QWEAS092- 15, QWEAS093-15, DFO097-11, RBCMI219- 14, KP976290, KP976289, KP976288, KP976287, KP976286, KJ125420, KJ125421, JQ267795 0‒2.39% Group 2 Remaining Mya arenaria (142 sequences) 11.55–13.94% 0‒1.99% Group 3 Mya truncata KF644116, KF644129 17.53–18.33% 15.54–17.53% 0.40% Group 4 Mya truncata KF644154, KF643769, KF643675, KF643403 17.53–18.73% 17.53–19.12% 7.17–11.55% 0 Group 5 Mya uzenensis KX534204 17.53–18.33% 16.33–17.93% 11.16–11.55% 13.15% ‒ Group 1 Group 2 Group 3 Group 4 Group 5 Group 1 All Mya japonica + Mya arenaria QWEAS092- 15, QWEAS093-15, DFO097-11, RBCMI219- 14, KP976290, KP976289, KP976288, KP976287, KP976286, KJ125420, KJ125421, JQ267795 0‒2.39% Group 2 Remaining Mya arenaria (142 sequences) 11.55–13.94% 0‒1.99% Group 3 Mya truncata KF644116, KF644129 17.53–18.33% 15.54–17.53% 0.40% Group 4 Mya truncata KF644154, KF643769, KF643675, KF643403 17.53–18.73% 17.53–19.12% 7.17–11.55% 0 Group 5 Mya uzenensis KX534204 17.53–18.33% 16.33–17.93% 11.16–11.55% 13.15% ‒ View Large Figure 5. View largeDownload slide Genetic distance range (A), phylogenetic trees (B, D) and candidate species inferred from Automatic Barcode Gap Discovery (ABGD) (C). A, the COI intraspecific distance values with M. truncata are 0, 0.33, 6.93 and 7.26% marked with green spots. B, phylogenetic tree inferred by Bayesian analysis (BI). C, ABGD, based on the COI alignment, with both models of Jukes-Cantor (JC69) and Kimura (K80) TS/TV. D, phylogenetic tree inferred by maximum likelihood (ML). Figure 5. View largeDownload slide Genetic distance range (A), phylogenetic trees (B, D) and candidate species inferred from Automatic Barcode Gap Discovery (ABGD) (C). A, the COI intraspecific distance values with M. truncata are 0, 0.33, 6.93 and 7.26% marked with green spots. B, phylogenetic tree inferred by Bayesian analysis (BI). C, ABGD, based on the COI alignment, with both models of Jukes-Cantor (JC69) and Kimura (K80) TS/TV. D, phylogenetic tree inferred by maximum likelihood (ML). Phylogenetic analyses For COI, 16S rDNA, 28S rDNA and the concatenated alignments, the TIM3 + I + G, TPM1uf + I, HKY and TVM + G evolutionary models were the best-fitted models, respectively. We report only the results of trees generated from COI sequences (Fig. 5B, D), which better resolved relationships at species level, and trees from the concatenated regions (Fig. 6). Trees constructed from 28S were uninformative with low support values at all levels and were omitted. Trees constructed from 16S were nearly identical in topology with the concatenated regions (except for differences in support values) and were also omitted. The resulting two consensus trees generated using BI and ML analyses were generally consistent; thus, a single topology was presented for the concatenated regions with support values indicated on branches (Fig. 6). This topology was in agreement with the COI trees and showed species-level groupings with high support (Fig. 6). Figure 6. View largeDownload slide Bayesian (BI) and maximum likelihood (ML) consensus tree inferred from concatenated sequences of COI, 16S and 28S genes. Figure 6. View largeDownload slide Bayesian (BI) and maximum likelihood (ML) consensus tree inferred from concatenated sequences of COI, 16S and 28S genes. In all COI trees, the ABGD Group 1 (all western Pacific soft-shell clams and two individuals of M. arenaria from the northeast Pacific) formed a clade with high node support (BI 1.00, ML 100%) and the ABGD Group 2 (remaining M. arenaria) grouped together with low to high node support before clustering to other species (BI 1.00, ML 42%, Fig. 5B, D). In the COI BI tree and the concatenated trees, M. arenaria clustered with M. japonica with high support, while M. truncata formed a sister clade with M. uzenensis with full node support (Figs 5B, 6). In the COI BI tree, M. japonica and M. arenaria QWEAS092-15 + QWEAS093-15 + DFO097-11 + RBCMI219-14 branched early within Mya with high node support (ML 95%), while the remaining M. arenaria clustered with the clade M. truncata + M. uzenensis with low support (ML 25%, Fig. 5D). In all COI trees, M. truncata were separated into two subclades, KF644116 + KF644129 and KF644154 + KF643769 + KF643675 + KF643403 (BI 0.86, ML 88). Divergence time estimates Assuming r1 = 0.35%, M. japonica would have split from M. arenaria ~9.4 Myr [95% highest posterior density (HPD) bounds: 6.7–12.5 Myr]. Mya uzenensis and M. truncata separated 8.1 Myr (95% HPD bounds: 4.0–11.9 Myr; Fig. 7). Time (T) to the most recent common ancestor for the Mya would be 9.9–20.0 Myr. Alternatively, assuming a faster substitution rate (r2 = 0.6% per lineage per Myr), M. japonica would have split from its congener, M. arenaria, ~5.7 Myr (95% HPD bounds: 4.1–7.5 Myr). Mya uzenensis and M. truncata separated 4.7 Mya (95% HPD bounds: 2.4–6.8 Myr). T for Mya would be 5.8–11.4 Myr (Fig. 7). Figure 7. View largeDownload slide Phylogenetic relationships among bivalve molluscs in Mya with estimates of divergence times, based on COI sequences. Bars in the tree correspond to 95% credibility intervals. Values above the bars were estimated by substitution rate r1 = 0.35% and below by substitution rate r2 = 0.6%. Scar bars on the bottom indicate time before present (Myr) derived from the two substitution rates (r1 and r2). Figure 7. View largeDownload slide Phylogenetic relationships among bivalve molluscs in Mya with estimates of divergence times, based on COI sequences. Bars in the tree correspond to 95% credibility intervals. Values above the bars were estimated by substitution rate r1 = 0.35% and below by substitution rate r2 = 0.6%. Scar bars on the bottom indicate time before present (Myr) derived from the two substitution rates (r1 and r2). DISCUSSION Mya japonica is a valid species Although the morphological differences between the clams identified as M. arenaria and M. japonica in the present study were minor, the ABGD delimitation, reciprocal monophyly, genetic divergence and spermatozoan ultramorphology of soft-shell clams from the western Pacific strongly support the recognition of M. japonica as a valid species. Within Mya, we found two clades, one comprising M. arenaria and M. japonica and the other containing both other Mya species examined. While M. arenaria and M. japonica were most related to each other, there was deep divergence between M. arenaria from the Northeast Pacific, Atlantic, Mediterranean and Barents Sea (Arctic) and M. japonica from the Pacific, which was further supported by ABGD. Furthermore, we found significant differences in spermatozoan ultramorphology between M. arenaria (Barents Sea) and M. japonica (Sea of Japan and Yellow Sea). The spermatozoan differences and genetic divergence reported in this study suggest long-term reproductive isolation between these two morphologically similar species. The fine morphology of spermatozoa from M. japonica examined in this study was very similar to M. japonica from Peter the Great Bay, Sea of Japan (Drozdov et al., 2009) and ‘M. arenaria oonogai’ from the coastal Korea Strait of Samcheonpo, Korea (Kim et al., 2011), which should all represent M. japonica s.s. Similarly, previous junior synonyms of M. arenaria described from the western Pacific, such as M. oonogai and M. arenaria oonogai, probably represent M. japonica. Shell morphological examinations and intraspecific variation of M. japonica Although individuals identified as M. arenaria and M. japonica were similar in appearance, there were recognizable differences in general shell shape, shell texture and appearance of the pallial line. Although the shape, size and orientation of the left valve chondrophore are often important diagnostic characters among Mya and related bivalves (Bernard, 1979), we found considerable variation in the characteristics of the chondrophore among all the material examined and it was not useful in discriminating between species. Given the comparative nature of the characteristics reported here (e.g. less/more impressed, elongated/shorter etc.) and overall degree of intraspecific variation, it may be difficult to identify individual clams using the descriptions contained herein, particularly for those who have not had the benefit of studying many specimens. Thus, we strongly recommend molecular or other analyses to confidently separate M. arenaria and M. japonica. Makiyama (1934, 1935) reported two forms of soft-shell clam in Japan, a northern form described as M. japonica and the southern form described as M. oonogai. The latter was referred to as M. japonica oonogai by Habe (1955) and Fujie (1957, 1962) and M. arenaria oonogai by Kwon et al. (2001), Kwon, Park & Lee (1993) and Yoo (1967). Fujie (1957) not only agreed that northern forms were M. japonica but also described a forma ‘α’, which was shorter and more ovate. Nagao & Inoue (1941) described differences between northern and southern Japanese forms and noted many intermediate specimens. The northern group, which possessed a thick, relatively short and nearly equilateral shell that was somewhat well rounded or slightly truncated posteriorly, was considered M. arenaria (due to synonymy with M. japonica). The southern group, which was characterized by its generally thin, long and somewhat inequilateral shell with a narrowly rounded or frequently acuminated posterior extremity, was considered an unnamed species. Nagao & Inoue (1941) also reported that the chondrophore in the northern group was nearly perpendicular to the antero-posterior diameter of the shell, and was narrower and more deeply excavated, with its inner border less convex than in the southern group. The posterior ridge was also not well differentiated from the chondrophore. Similarly, Scarlato (1981) noted two morphological forms of M. japonica in Russian Far Eastern seas. One form is distributed in low-boreal (temperate) waters (Aniva Bay in southern Sakhalin, southern Kurile Islands and Peter the Great Bay), where this species reaches its maximum size and possesses the typical shape with a smooth surface and fine, regular growth lines. The other form was found in the northern part of its distributional range (Bering and Chukotsk seas, eastern Kamchatka and northern Sea of Okhotsk), where M. japonica never reaches its maximum length and often has an atypical shape, sometimes appearing misshapen or deformed. However, our examination of M. japonica from Peter the Great Bay indicates that these variations also occur in the same geographical locality and are related to different environments (Fig. 2) along with intermediate specimens. The co-occurrence of different morphological forms and intermediate specimens in a single geographical region strongly suggests that these forms are ecophenotypes. It should be noted that the ‘typical’ thin and regular form reported here fully complies with ‘M. japonica oonogai’ collected from Hokkaido and Nagasaki, Japan, and our ‘thick-shelled’ form is very similar to M. japonica forma ‘α’, as illustrated by Fujie (1957). In addition, the M. japonica of Fujie (1957) is intermediate between the two forms we reported here. The results of our phylogenetic and morphological analyses indicate that M. japonica is a single species, which displays a high degree of phenotypic plasticity in shell morphology, and many previously reported forms/species combinations probably reflect ecophenotypic variation or represent natural intraspecific variation. Distribution of M. arenaria and M. japonica The soft-shell clam M. arenaria is widely accepted to have a circumboreal distribution (Bernard, 1979, 1983; Strasser, 1999; Coan et al., 2000; Huber, 2010; Zhang et al., 2012). Previous molecular studies indicated a complex colonization history between Europe and North America, but these analyses lacked samples from the northwest Pacific – the putative range of M. japonica (Strasser & Barber, 2009; Layton et al., 2014; Barco et al., 2016; Cross et al., 2016; Lasota et al., 2016). Our genetic and spermatozoan ultramorphological analyses strongly suggest that soft-shell clams in the northwest Pacific represent M. japonica. Two publicly available sequences identified as ‘M. arenaria’ (the valid species name at the time of identification) from China (QWEAS092-15 + QWEAS093-15) were highly similar (< 2% difference COI; Supporting Information, Table S2) to western Pacific M. japonica. These individuals are probably M. japonica and do not reflect the presence of M. arenaria in the western Pacific. MacNeil (1965) reported that M. arenaria occurred in Japan only in the middle Miocene and that later and recent specimens from the Okhotsk Sea, west coast of Japan and China are M. japonica, which was supported by Lutaenko & Noseworthy (2012) and Scarlato (1981). Numerous other authors have also reported M. japonica from China, Japan, Korea and the Russian Far Eastern seas (Jay, 1856; Makiyama, 1935; Habe, 1955, 1977; Tchang, Qi & Li, 1955; Fujie, 1957; Scarlato, 1981; Okutani, 2000). While both M. arenaria and M. japonica were present along western North America at one time, there is no palaeontological or archaeological evidence that either species survived the Pleistocene glaciation, with the possible exception of relict populations in the Bering Sea (Carlton, 1979). Mya arenaria was first introduced into the northeastern Pacific via Eastern Oyster Crassostrea virginica Gmelin, 1791 seed plantings in San Francisco Bay, California, during the 1860–1870s (Carlton, 1979). Since that time, it has spread northward, by either intentional plantings or natural dispersal, to southern Alaska (Carlton, 1979; Strasser, 1999; Powers et al., 2006). In contrast, until the present study, there were no recent eastern Pacific records of M. japonica or any evidence of intentional or accidental introductions of western Pacific Mya into the eastern Pacific. Two publicly available sequences identified as ‘M. arenaria’ from the British Columbia, one collected from Quascilla Bay (DFO097-11) and the other from Graham Island (RBCMI219-14), were highly similar (< 2% difference COI; Supporting Information, Table S2) to all western Pacific M. japonica. These individuals, both collected in June 2011, are probably M. japonica and went unnoticed because of morphological similarities to M. arenaria, the expected species in British Columbia. Thus, these records represent the first modern documented occurrence of M. japonica in the northeast Pacific. In addition, after examining photographs of these two clams, the authors are confident in classifying them as M. japonica based on the morphological criteria used in this study. We are uncertain of the source of introduction or status of possible eastern Pacific M. japonica populations. However, similar to M. arenaria, M. japonica was probably introduced via oyster plantings, but in this case, the Pacific Oyster Crassostrea gigas Thunberg, 1793. These plantings, which occurred from the early 1910s until after the Second World War (Carlton, 1979; Lavoie, 2005), are also thought to have led to other bivalve introductions, such as the Japanese Littleneck Ruditapes philippinarum Adams & Reeve, 1850 and Quadrate Trapezium Neotrapezium liratum Reeve, 1843. Alternatively, M. japonica may have been introduced by ballast water, as is speculated with the Asian Semele Theora lubrica Gould, 1861 and Asian brackish-water clam P. amurensis in California (Fofonoff et al., 2003). In either case, M. japonica was probably introduced decades after M. arenaria became widespread and well established in the eastern Pacific (Carlton, 1979), which reduced the likelihood of detection of M. japonica in areas where M. arenaria was already present or in nearby localities. Furthermore, both individuals were larger-sized adults collected nearly 400 km apart, which suggests there may be viable populations of M. japonica in at least British Columbia or possibly other Pacific Northeast locations. Albeit unlikely, the presence of individuals with high COI similarly to M. japonica could represent F1 hybrids with maternal M. japonica contribution or possibly reflect historical or cryptic hybridization (Pfenninger, Reinhardt & Streit, 2002; Rees, Dioli & Kirkendall, 2003) between M. arenaria and M. japonica. While there were significant differences in sperm morphology between M. arenaria and M. japonica, we are uncertain if these differences (or ecological/behavioural differences) would prohibit hybridization. According to Laursen (1966), historical reports of M. arenaria from the Arctic Ocean are erroneous and should represent M. truncata ovata Jensen, 1900, which is currently recognized as a junior synonym of M. pseudoarenaria (see Huber, 2010). However, we identified M. arenaria from the Russian waters of the Barents Sea, a marginal sea of the Arctic Ocean, and it has also been reported further west near Forsøl, Norway (Crocetta & Turolla, 2011). A rigorous genetic analysis (incorporating both nuclear and mitochondrial markers) of M. arenaria and M. japonica collected from across the eastern and western Pacific, Bering Sea, Arctic Ocean and greater Atlantic is required to fully resolve the zoogeography of both species, determine the extent of M. japonica introductions in the eastern Pacific, identify any possible cryptic introductions of M. arenaria in the western Pacific, determine the identity of soft-shell clams in the Bering Sea and examine possible hybridization of M. arenaria and M. japonica. Divergence time and evolutionary/migration history We estimate that M. japonica diverged from M. arenaria 4.1–12.5 Myr; this agrees with the evolution of M. arenaria, which evolved during early Pliocene to late Miocene from its ancestor Mya fujiei MacNeil, 1965, which arose during middle Miocene (MacNeil, 1965; Strauch, 1972). Mya arenaria first appeared in the fossil record from Japan in late Miocene formations and almost at the same time from the eastern Pacific in California (Strauch, 1972). Mya japonica was also probably present in the eastern Pacific at least during the Pleistocene (MacNeil, 1965; Carlton, 1979). Although the formation of the Bering Strait during the Pliocene overlapped with the divergence of M. arenaria and M. japonica, it is likely that only M. arenaria crossed the Arctic to reach the western and eastern Atlantic (MacNeil, 1965; Bernard, 1979; Vermeij, 1989). Strauch (1972) suggested M. arenaria may also have migrated through the Central American Passage from the Pacific to the Atlantic, where it then spread to Europe in the late Pliocene (Strasser, 1999), but fossil evidence for this migratory route is lacking (Bernard, 1979). During the Pleistocene glaciation, M. arenaria was extirpated from all areas except the western Atlantic and M. japonica was extirpated from all areas except the western Pacific (Strasser, 1999). However, the status of fossil, archaeological and recent soft-shell clams in the Bering Sea requires further investigation, which could represent relict populations of M. arenaria or M. japonica, adventitious individuals of M. japonica from the western Pacific or an undescribed cryptic species (MacNeil, 1965; Bernard, 1979; Carlton, 1979). The evolutionary history of Mya is complex and there is much uncertainty regarding the identification and evolutionary relationships of fossil material (MacNeil, 1965) and their connection with extant forms, particularly as it relates to the synonymy of extant species and species identified solely from fossil material (Petersen, 1999). Due to the morphological similarities between M. arenaria and M. japonica and levels of phenotypic plasticity, it may be difficult to fully resolve their evolutionary and migration history. MacNeil (1965) and Strauch (1972) reported that the oldest records of M. truncata and Mya cuneiformis Böhm 1915 were from the middle Miocene, both of which were direct descendants of the late Oligocene or early Miocene species, Mya salmonensis Clark, 1932, while M. pseudoarenaria and Mya priapus Tilesius, 1822 derived from M. cuneiformis in the late Miocene. Nakashima (1999) determined that the oldest records of M. pseudoarenaria, M. cuneiformis and M. truncata were from the early Miocene and M. uzenensis may have directly evolved from M. salmonensis. Our results indicate that the separation time of M. uzenensis and M. truncata (3.9–21.0 Myr) is nearly the same as M. japonica and M. arenaria and that M. uzenensis and M. truncata share a common ancestor, but it is unclear if M. uzenensis directly evolved from M. salmonensis or M. cuneiformis. Mya truncata may reflect a species complex The ABGD analysis indicated that M. truncata was a candidate species complex (Fig. 5C), which included two distinct groups (ABGD Groups 3 and 4; Table 3). The COI genetic distances between individuals of these groups ranged from 6.93 to 7.26% (Table 3). While these are intermediate values between the maximum intraspecific p-distance (1.70%) within M. japonica + M. arenaria and the minimum interspecific p-distance (10.68%) among the other Mya species, this may reflect a true species-level difference, which is consistent with Petersen (1999) and Layton et al. (2014). At present, we are unable to fully resolve the taxonomic status of M. truncata. All the M. truncata COI sequences examined in this study were obtained from individuals collected in Nunavut, Canada (Layton et al., 2014). However, individuals in ABGD Group 3 were collected near Cornwallis Island, whereas individuals in ABGD Group 4 were collected around Igloolik Island, indicating possible spatial structure to the putative species complex. Specimens collected from the entire range of M. truncata (eastern Pacific, Arctic and northern Atlantic) and additional data (e.g. spermatozoan ultramorphology) are needed to clarify the taxonomic status of M. truncata. CONCLUSIONS Our results indicate that (1) M. japonica is a valid species, distributed in the northern Pacific, with introduced populations of unknown size and distribution in the northeast Pacific; (2) M. arenaria is distributed in the northeast Pacific, northern Atlantic, Mediterranean and Barents Sea (Arctic); (3) there were significant species-level differences in the spermatozoan ultramorphology of M. arenaria and M. japonica; (4) differences in shell morphology were minor, and given the levels of intraspecific variation noted among varying habitats, shell morphology alone may not be an effective method to differentiate M. arenaria and M. japonica; (5) the congeners, M. arenaria, M. japonica, M. uzenensis and M. truncata were all differentiated by DNA barcoding, with evidence of cryptic diversity in M. truncata from the Canadian Arctic; and (6) the timing of divergence between M. arenaria and M. japonica supports the evolutionary hypothesis that M. arenaria evolved during early Pliocene to late Miocene. Further collections of Mya species from across their entire range (especially Asia and Arctic Ocean/Bering Sea where collections have been scarce) and a rigorous analysis incorporating an integrative taxonomic approach using ecological, molecular, morphological and/or other biological data are required to fully resolve taxonomic relationships and distributions of these highly morpho-variable species. SUPPORTING INFORMATION Supporting information may be found in the online version of this article at the publisher’s web-site: Table S1. List of all specimens used for molecular analyses. Table S2. Interspecific and intraspecific uncorrected pairwise distances at COI among species of Mya and Potamocorbula. Table S3. Interspecific and intraspecific uncorrected pairwise distances at 16S among species of Mya and Potamocorbula. ACKNOWLEDGEMENTS This study was supported by the Special Funds for the Young Scholars of Taxonomy of the Chinese Academy of Sciences (No. ZSBR-009), the Bureau of International Cooperation Chinese Academy of Sciences, the Far Eastern Branch of the Russian Academy of Sciences (No. 15-I-6–007 o) and the Russian Foundation for Basic Research (No. RFFI-ofi 15-29-0245615). The authors are grateful to Amy Driskell (Smithsonian Laboratories of Analytical Biology) for sequencing the Chesapeake Bay Mya arenaria samples. We thank Prof. Jim Carlton (Williams College, USA) for helpful suggestions on the manuscript. Irina E. Volvenko (Zoological Museum, Far Eastern Federal University, Russia) kindly took photographs of shells. REFERENCES Aguilar R , Ogburn MB , Driskell AC , Weigt LA , Groves MC , Hines AH . 2016 . Gutsy genetics: identification of digested piscine prey items in the stomach contents of sympatric native and introduced warmwater catfishes via DNA barcoding . Environmental Biology of Fishes 100 : 325 – 336 . Google Scholar CrossRef Search ADS Alfaro ME , Zoller S , Lutzoni F . 2003 . Bayes or bootstrap? A simulation study comparing the performance of Bayesian Markov chain Monte Carlo sampling and bootstrapping in assessing phylogenetic confidence . Molecular Biology and Evolution 20 : 255 – 266 . Google Scholar CrossRef Search ADS PubMed Barco A , Raupach MJ , Laakmann S , Neumann H , Knebelsberger T . 2016 . Identification of North Sea molluscs with DNA barcoding . Molecular Ecology Resources 16 : 288 – 297 . Google Scholar CrossRef Search ADS PubMed Bernard FR . 1979 . Identification of the living Mya (Bivalvia: Myoida) . Venus: Japanese Journal of Malacology 38 : 185 – 204 . Bernard FR . 1983 . Catalogue of the living Bivalvia of the eastern Pacific Ocean: Bering Strait to Cape Horn . Ottawa : Love Printing Services Ltd . Bouchet P , Gofas S . 2013 . Mya Linnaeus, 1758 . In: MolluscaBase (2017) . World Register of Marine Species . Available at: http://marinespecies.org/aphia.php?p=taxdetails&id=138211 (accessed 10 May 2017 ). Bouckaert R , Heled J , Kühnert D , Vaughan T , Wu CH , Xie D , Suchard MA , Rambaut A , Drummond AJ . 2014 . BEAST 2: a software platform for Bayesian evolutionary analysis . PLoS Computational Biology 10 : e1003537 . Google Scholar CrossRef Search ADS PubMed Carlton JT . 1979 . History, biogeography, and ecology of the introduced marine and estuarine invertebrates of the Pacific coast of North America . Unpublished Ph.D. Thesis, University of California . Carmona L , Lei BR , Pola M , Gosliner TM , Valdés Á , Cervera JL . 2014 . Untangling the Spurilla neapolitana (Delle Chiaje, 1841) species complex: a review of the genus Spurilla Bergh, 1864 (Mollusca: Nudibranchia: Aeolidiidae) . Zoological Journal of the Linnean Society 170 : 132 – 154 . Google Scholar CrossRef Search ADS Carstensen D , Laudien J , Leese F , Arntz W , Held C . 2009 . Genetic variability, shell and sperm morphology suggest that the surf clams Donax marincovichi and D. obesulus are one species . Journal of Molluscan Studies 75 : 381 – 390 . Google Scholar CrossRef Search ADS Coan EV , Scott PV , Bernard FR . 2000 . Bivalve seashells of western North America: marine bivalve mollusks from Arctic Alaska to Baja California. Santa Barbara Museum of Natural History Monographs Number 2; Studies in Biodiversity Number 2 . Santa Barbara : Santa Barbara Museum of Natural History . Congleton WR , Vassiliev T , Bayer RC , Pearce BR , Jacques J , Gillman C . 2006 . Trends in Maine softshell clam landings . Journal of Shellfish Research 25 : 475 – 480 . Google Scholar CrossRef Search ADS Crocetta F , Turolla E . 2011 . Mya arenaria Linné, 1758 (Mollusca: Bivalvia) in the Mediterranean Sea: its distribution revisited . Journal of Biological Research-Thessaloniki 16 : 188 – 193 . Cross ME , Bradley CR , Cross TF , Culloty S , Lynch S , McGinnity P , O’Riordan RM , Vartia S , Prodohl PA . 2016 . Genetic evidence supports recolonisation by Mya arenaria of western Europe from North America . Marine Ecology Progress Series 549 : 99 – 112 . Google Scholar CrossRef Search ADS Darriba D , Taboada GL , Doallo R , Posada D . 2012 . jModelTest 2: more models, new heuristics and parallel computing . Nature Methods 9 : 772 . Google Scholar CrossRef Search ADS PubMed Drozdov AL , Sharina SN , Tyurin SA . 2009 . Sperm ultrastructure in representatives of six bivalve families from Peter the Great Bay, Sea of Japan . Russian Journal of Marine Biology 35 : 236 – 241 . Google Scholar CrossRef Search ADS Dzyuba SM , Maslennikova LA . 1987 . Reproductive-cycle of the bivalve mollusk Mya japonica in Peter the Great Bay of the Sea of Japan . Biologiya Morya-Marine Biology 2 : 38 – 41 . Eckelbarger KJ , Young CM . 1999 . Ultrastructure of gametogenesis in a chemosynthetic mytilid bivalve (Bathymodiolus childressi) from a bathyal, methane seep environment (northern Gulf of Mexico) . Marine Biology 135 : 635 – 646 . Google Scholar CrossRef Search ADS Edgecombe GD , Giribet G , Wheeler WC . 2002 . Phylogeny of Henicopidae (Chilopoda: Lithobiomorpha): a combined analysis of morphology and five molecular loci . Systematic Entomology 27 : 31 – 64 . Google Scholar CrossRef Search ADS Elderkin CL , Clewing C , Ndeo OW , Albrecht C . 2016 . Molecular phylogeny and DNA barcoding confirm cryptic species in the African freshwater oyster Etheria elliptica Lamarck, 1807 (Bivalvia: Etheriidae) . Biological Journal of the Linnean Society 118 : 369 – 381 . Google Scholar CrossRef Search ADS Fofonoff PW , Ruiz GM , Steves B , Carlton JT . 2003 . California Non-native Estuarine and Marine Organisms (Cal-NEMO) system . Available at: http://invasions.si.edu/nemesis (accessed 2 May 2017 ). Folmer O , Black M , Hoeh W , Lutz R , Vrijenhoek R . 1994 . DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates . Molecular Marine Biology and Biotechnology 3 : 294 – 299 . Google Scholar PubMed Fujie T . 1957 . On the myarian pelecypods of Japan: part I. Summary of the study of the genus Mya from Hokkaido . Journal of the Faculty of Science, Hokkaido University. Series 4, Geology and Mineralogy 9 : 381 – 413 . Fujie T . 1962 . On the myarian pelecypods of Japan: Part II. Geological and geographical distribution of fossil and recent species, genus Mya . Journal of the Faculty of Science, Hokkaido University. Series 4, Geology and Mineralogy 11 : 399 – 429 . Geraghty J , Russell MP , Dollahon N . 2008 . A quantitative assessment of spermatozoan morphology in Nutricola confusa and Nutricola tantilla (Bivalvia: Veneridae) . Veliger 50 : 263 – 268 . Gomoiu M-T , Alexandrov B , Shadrin N , Zaitsev Y . 2002 . The Black Sea – a recipient, donor and transit area for alien species . In: Leppäkoski E , Gollasch S , Olenin S , eds. Invasive aquatic species of Europe. Distribution, impacts and management . Dordrecht : Kluwer Academic Publishers , 341 – 350 . Google Scholar CrossRef Search ADS Guindon S , Dufayard JF , Lefort V , Anisimova M , Hordijk W , Gascuel O . 2010 . New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0 . Systematic Biology 59 : 307 – 321 . Google Scholar CrossRef Search ADS PubMed Gwo JC , Liou CH , Cheng CH . 1996 . Ultrastructure of the spermatozoa of the Pacific oyster Crassostrea gigas (Mollusca, Bivalvia, Ostreidae) . Journal of Submicroscopic Cytology and Pathology 28 : 395 – 400 . Habe T . 1955 . Fauna of Akkeshi Bay. XXI. Pelecypoda and Scaphopoda . Publications of the Akkeshi Marine Biological Station 4 : 1 – 31 . Habe T . 1977 . Systematics of Mollusca in Japan. Bivalvia and Scaphopoda . Tokyo : Hokuryukan . Hillis DM , Bull JJ . 1993 . An empirical test of bootstrapping as a method for assessing confidence in phylogenetic analysis . Systematic Biology 42 : 182 – 192 . Google Scholar CrossRef Search ADS Hodgson AN , Bernard RTF , Vanderhorst G . 1990 . Comparative spermatology of three species of Donax (Bivalvia) from South-Africa . Journal of Molluscan Studies 56 : 257 – 265 . Google Scholar CrossRef Search ADS Hodgson AN , Devilliers CJ , Bernard RTF . 1987 . Comparative spermatology of two morphologically similar species of Solen (Mollusca, Bivalvia) . South African Journal of Zoology 22 : 264 – 268 . Google Scholar CrossRef Search ADS Huber M . 2010 . Compendium of bivalves. A full-color guide to 3’300 of the world’s marine bivalves: a status on Bivalvia after 250 years of research . Hackenheim : ConchBooks . Introini GO , Passos FD , Recco-Pimentel SM . 2013 . Comparative study of sperm ultrastructure of Donax hanleyanus and Donax gemmula (Bivalvia: Donacidae) . Acta Zoologica 94 : 261 – 266 . Google Scholar CrossRef Search ADS ITIS . 2017 . Mya Linnaeus, 1758 . Integrated Taxonomic Information System . Available at: https://www.itis.gov/servlet/SingleRpt/SingleRpt?search_topic=TSN&search_value=81691#null (accessed 10 May 2017 ). Jay JC . 1856 . Report on the shells collected by the Japan Expedition together with a list of Japan shells . In: Perry MC , Hawks FL , eds. Narrative of the expedition of an American squadron to the China Sea and Japan, performed in the year 1852, 1853 and 1854, under the command of Commodore M. C. Perry, United States Navy, by order of the Government of the United States . Washington : Tucker , 289 – 297 . Katoh K , Standley DM . 2013 . MAFFT multiple sequence alignment software version 7: improvements in performance and usability . Molecular Biology and Evolution 30 : 772 – 780 . Google Scholar CrossRef Search ADS PubMed Kienberger K , Carmona L , Pola M , Padula V , Gosliner TM , Cervera JL . 2016 . Aeolidia papillosa (Linnaeus, 1761) (Mollusca: Heterobranchia: Nudibranchia), single species or a cryptic species complex? a morphological and molecular study . Zoological Journal of the Linnean Society 177 : 481 – 506 . Google Scholar CrossRef Search ADS Kim JH , Chung JS , Park YJ . 2011 . Ultrastructures of germ cells during spermatogenesis and taxonomic values in sperm morphology in male Mya arenaria oonogai (Heterodonta: Myidae) . Korean Journal of Malacology 27 : 377 – 386 . Google Scholar CrossRef Search ADS Kwon OK , Min DK , Lee JR , Lee JS , Je JG , Choe BL . 2001 . Korean mollusks with color illustration . Seoul : Hanguel Publishing Company . Kwon OK , Park GM , Lee JS . 1993 . Colored shells of Korea . Seoul : Academy Publishing Co . Lasota R , Pierscieniak K , Garcia P , Simon-Bouhet B , Wolowicz M . 2016 . Large-scale mitochondrial COI gene sequence variability reflects the complex colonization history of the invasive soft-shell clam, Mya arenaria (L.) (Bivalvia) . Estuarine, Coastal and Shelf Science 181 : 256 – 265 . Google Scholar CrossRef Search ADS Laursen D . 1966 . The genus Mya in the Arctic Region . Malacologia 3 : 399 – 418 . Lavoie R . 2005 . Oyster culture in North America history, present and future . The 1st International Oyster Symposium Proceedings 19 : 9 . Layton KK , Martel AL , Hebert PD . 2014 . Patterns of DNA barcode variation in Canadian marine molluscs . PLoS ONE 9 : e95003 . Google Scholar CrossRef Search ADS PubMed Lutaenko KA , Noseworthy RG . 2012 . Catalogue of the living Bivalvia of the continental coast of the Sea of Japan (East Sea) . Vladivostok : Dalnauka . MacNeil FS . 1965 . Evolution and distribution of the genus Mya, and Tertiary migrations of Mollusca . United States, Department of the Interior, Geological Survey, Professional Paper 483G : 1 – 51 . Makiyama J . 1934 . The Asagaian molluscs of Yotukara and Matchgar . Memoirs of the College of Science Kyoto Imperial University, Series B 10 : 156 – 159 . Makiyama J . 1935 . The fossils of the genus Mya . Warera no Kobutsu 4 : 135 – 139 . Marko PB . 2002 . Fossil calibration of molecular clocks and the divergence times of geminate species pairs separated by the Isthmus of Panama . Molecular Biology and Evolution 19 : 2005 – 2021 . Google Scholar CrossRef Search ADS PubMed Nagao T , Inoue T . 1941 . Myarian fossils from the Cenozoic deposits of hokkaidô and Karahuto . Journal of the Faculty of Science, Hokkaido Imperial University. Series 4, Geology and Mineralogy 6 : 143 – 158 . Nakashima R . 1999 . Mya (Bivalvia: Myidae) from the Upper Miocene and Pliocene formations in Hokkaido, Northern Japan . Venus: Japanese Journal of Malacology 58 : 201 – 216 . Navarro D . 2014 . Learning statistics with R: a tutorial for psychology students and other beginners, Version 0.5 . Adelaide : University of Adelaide . Nikula R , Strelkov P , Väinölä R . 2007 . Diversity and trans-arctic invasion history of mitochondrial lineages in the North Atlantic Macoma balthica complex (Bivalvia: Tellinidae) . Evolution 61 : 928 – 941 . Google Scholar CrossRef Search ADS PubMed Okutani T . 2000 . Marine Mollusks in Japan . Tokyo : Tokai University Press . Padial JM , Miralles A , De la Riva I , Vences M . 2010 . The integrative future of taxonomy . Frontiers in Zoology 7 : 16 . Google Scholar CrossRef Search ADS PubMed Petersen GH . 1999 . Five recent Mya species, including three new species and their fossil connections . Polar Biology 22 : 322 – 328 . Google Scholar CrossRef Search ADS Pfenninger M , Reinhardt F , Streit B . 2002 . Evidence for cryptic hybridization between different evolutionary lineages of the invasive clam genus Corbicula (Veneroida, Bivalvia) . Journal of Evolutionary Biology 15 : 818 – 829 . Google Scholar CrossRef Search ADS Porter RG . 1974 . Reproductive cycle of the soft-shell clam, Mya arenaria, at Skagit Bay, Washington . Fishery Bulletin 72 : 648 – 656 . Powers SP , Bishop MA , Grabowski JH , Peterson CH . 2006 . Distribution of the invasive bivalve Mya arenaria L. on intertidal flats of southcentral Alaska . Journal of Sea Research 55 : 207 – 216 . Google Scholar CrossRef Search ADS Puillandre N , Lambert A , Brouillet S , Achaz G . 2012 . ABGD, Automatic Barcode Gap Discovery for primary species delimitation . Molecular Ecology 21 : 1864 – 1877 . Google Scholar CrossRef Search ADS PubMed R Core Team . 2015 . R: a language and environment for statistical computing . Vienna : R Foundation for Statistical Computing . Ratnasingham S , Hebert PD . 2007 . Bold: the barcode of life data system (http://www.barcodinglife.org) . Molecular Ecology Notes 7 : 355 – 364 . Google Scholar CrossRef Search ADS PubMed Rees DJ , Dioli M , Kirkendall LR . 2003 . Molecules and morphology: evidence for cryptic hybridization in African Hyalomma (Acari: Ixodidae) . Molecular Phylogenetics and Evolution 27 : 131 – 142 . Google Scholar CrossRef Search ADS PubMed Reynolds ES . 1963 . The use of lead citrate at high pH as an electron-opaque stain in electron microscopy . The Journal of Cell Biology 17 : 208 – 212 . Google Scholar CrossRef Search ADS PubMed Ronquist F , Teslenko M , van der Mark P , Ayres DL , Darling A , Höhna S , Larget B , Liu L , Suchard MA , Huelsenbeck JP . 2012 . MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space . Systematic Biology 61 : 539 – 542 . Google Scholar CrossRef Search ADS PubMed Scarlato OA . 1981 . Bivalve mollusks of temperate latitudes of the western portion of the Pacific Ocean. Guide books on the fauna of the USSR . Leningrad : Zoological Institute, USSR Academy of Sciences . Schlick-Steiner BC , Steiner FM , Seifert B , Stauffer C , Christian E , Crozier RH . 2010 . Integrative taxonomy: a multisource approach to exploring biodiversity . Annual Review of Entomology 55 : 421 – 438 . Google Scholar CrossRef Search ADS PubMed Schneider CA , Rasband WS , Eliceiri KW . 2012 . NIH Image to ImageJ: 25 years of image analysis . Nature Methods 9 : 671 – 675 . Google Scholar CrossRef Search ADS PubMed Strasser M . 1999 . Mya arenaria – an ancient invader of the North Sea coast . Helgoländer Meeresuntersuchungen 52 : 309 – 324 . Google Scholar CrossRef Search ADS Strasser CA , Barber PH . 2009 . Limited genetic variation and structure in softshell clams (Mya arenaria) across their native and introduced range . Conservation Genetics 10 : 803 – 814 . Google Scholar CrossRef Search ADS Strauch F . 1972 . Phylogenese, adaptation und migration einiger nordischer mariner Molluskengenera (Neptunea, Panomya, Cyrtodaria und Mya) . Abhandlungen der Senckenbergischen Naturforschenden Gesellschaft 531 : 1 – 211 . Streftaris N , Zenetos A . 2006 . Alien marine species in the Mediterranean – the 100 ‘worst invasives’ and their impact . Mediterranean Marine Science 7 : 87 – 118 . Google Scholar CrossRef Search ADS Tamura K , Stecher G , Peterson D , Filipski A , Kumar S . 2013 . MEGA6: molecular evolutionary genetics analysis version 6.0 . Molecular Biology and Evolution 30 : 2725 – 2729 . Google Scholar CrossRef Search ADS PubMed Tchang S , Qi Z , Li J . 1955 . Economic marine Mollusca in North China Sea . Beijing : Science Press . Torroglosa ME , Gimenez J . 2015 . Sperm ultrastructure in two species of Brachidontes (Bivalvia, Mytilidae) from the south-western Atlantic Ocean . Journal of the Marine Biological Association of the United Kingdom 95 : 991 – 998 . Google Scholar CrossRef Search ADS Vaidya G , Lohman DJ , Meier R . 2011 . SequenceMatrix: concatenation software for the fast assembly of multi-gene datasets with character set and codon information . Cladistics 27 : 171 – 180 . Google Scholar CrossRef Search ADS Väinölä R . 2003 . Repeated trans-Arctic invasions in littoral bivalves: molecular zoogeography of the Macoma balthica complex . Marine Biology 143 : 935 – 946 . Google Scholar CrossRef Search ADS Vermeij GJ . 1989 . Invasion and extinction: the last three million years of North Sea pelecypod history . Conservation Biology 3 : 274 – 281 . Google Scholar CrossRef Search ADS Whiting MF , Carpenter JC , Wheeler QD , Wheeler WC . 1997 . The Strepsiptera problem: phylogeny of the holometabolous insect orders inferred from 18S and 28S ribosomal DNA sequences and morphology . Systematic Biology 46 : 1 – 68 . Google Scholar PubMed Xiong B , Kocher TD . 1991 . Comparison of mitochondrial DNA sequences of seven morphospecies of black flies (Diptera: Simuliidae) . Genome 34 : 306 – 311 . Google Scholar CrossRef Search ADS PubMed Yoo JS . 1967 . Korean shells in colour . Seoul : Il Ji Sa Publishing Co . Zhang J , Xu F , Liu R . 2012 . The Myidae (Mollusca, Bivalvia) from Chinese waters with description of a new species . Zootaxa 3383 : 39 – 60 . Google Scholar CrossRef Search ADS © 2018 The Linnean Society of London, Zoological Journal of the Linnean Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Zoological Journal of the Linnean Society Oxford University Press

A tale of two soft-shell clams: an integrative taxonomic analysis confirms Mya japonica as a valid species distinct from Mya arenaria (Bivalvia: Myidae)

Loading next page...
 
/lp/ou_press/a-tale-of-two-soft-shell-clams-an-integrative-taxonomic-analysis-v7ywrzjSWy
Publisher
The Linnean Society of London
Copyright
© 2018 The Linnean Society of London, Zoological Journal of the Linnean Society
ISSN
0024-4082
eISSN
1096-3642
D.O.I.
10.1093/zoolinnean/zlx107
Publisher site
See Article on Publisher Site

Abstract

Abstract The soft-shell clam Mya arenaria Linnaeus, 1758 is a commercially important fishery resource that occurs in boreal and temperate environments in the Northern Hemisphere. Whether the soft-shell clam is a single species with a circumboreal range or a species complex also comprising Mya japonica Jay, 1857 distributed in the north Pacific has long been debated by malacologists and palaeontologists based on slight differences in shell morphology. We used an integrative taxonomic approach incorporating available Mya spp. mitochondrial COI and 16S rRNA, and nuclear 28S rRNA gene sequences, as well as spermatozoan and shell morphological characters to test the validity of M. japonica and examine the range of soft-shell clam distribution. Although differences in shell morphology were minor, the results from tree topologies, pairwise uncorrected p-distances, Automatic Barcode Gap Discovery (ABGD) and spermatozoan ultramorphological data confirm the validity of M. japonica in both its endemic region in the northwest Pacific, and as here newly reported introduced populations in British Columbia in the northeast Pacific, and show that M. arenaria is distributed in the northeast Pacific, North Atlantic, Barents Sea (Arctic Ocean) and Mediterranean. We estimate these two closely related sister species diverged 4.1–12.5 Myr during early Pliocene to late Miocene, which is consistent with current evolutionary theory regarding M. arenaria. In addition, ABGD indicated the congener Mya truncata Linnaeus, 1758 may represent a species complex, but additional evidence is still needed to clarify its taxonomic status. biogeography, cryptic species, DNA barcoding, molecular phylogeny, Mollusca, spermatozoan ultramorphology INTRODUCTION The soft-shell clam Mya arenaria Linnaeus, 1758, described from Europe, is an edible bivalve currently found in boreal and northern temperate marine and estuarine environments with a complex, but widely accepted migration history (MacNeil, 1965; Laursen, 1966; Strauch, 1972; Strasser, 1999). Mya arenaria originated in the western Pacific probably around Japan during the Miocene, spreading to the eastern Pacific later in the Miocene and then to both coasts of the Atlantic Ocean during the Pliocene (MacNeil, 1965). It was extirpated in the eastern Pacific and eastern Atlantic during the Pleistocene glaciation, surviving only in the western Atlantic with uncertain status in the Bering Sea and western Pacific (MacNeil, 1965; Carlton, 1979; Strasser, 1999). It re-invaded the eastern Atlantic during the Middle Ages and the eastern Pacific during the late 1800s, both invasions mediated by anthropogenic activities (Fujie, 1957, 1962; MacNeil, 1965; Bernard, 1979; Strasser, 1999; Cross et al., 2016). It has more recently appeared in the Mediterranean and Black Sea (Strasser, 1999; Gomoiu et al., 2002). Historical Arctic records are now generally regarded as belonging to a congener, Mya pseudoarenaria Schlesch, 1931 (Laursen, 1966). Mya arenaria is a major fishery resource in the western Atlantic (Congleton et al., 2006). However, in many parts of its introduced range M. arenaria can negatively impact native benthic communities (Streftaris & Zenetos, 2006; Crocetta & Turolla, 2011), disrupt native fisheries (Carlton, 1979) and represent a significant portion of benthic biomass (Strasser, 1999). Within the genus Mya there is much disagreement regarding evolutionary history, species distributions and taxonomy (MacNeil, 1965; Bernard, 1979; Petersen, 1999; Huber, 2010; Bouchet & Gofas, 2013; ITIS, 2017), particularly regarding the number of extant species and the relationships between extant forms and species described from fossil materials (Petersen, 1999). The most extensively studied species (both ecologically and taxonomically) is M. arenaria, but questions regarding taxonomic status and geographical range remain. Mya arenaria is reported to have numerous synonyms: M. communis Megerle von Mühlfeld, 1811; M. lata Sowerby, 1815; M. acuta Say, 1822; M. acuta mercenaria Say, 1822; M. subovata Woodward, 1833; M. subtruncata Woodward, 1833; M. alba Agassiz, 1839; M. corpulenta Conrad, 1845; M. japonica Jay, 1856; Sphenia ovoidea Carpenter, 1864; M. hemphilli Newcomb, 1874; M. declivis Pennant, 1777; M. elongata Locard, 1886; M. arenaria var. ovata Jensen, 1900; M. paternalis Matsumoto, 1930; M. arenaria oonogai Makiyama, 1935; M. japonica oonogai Makiyama, 1935; M. oonogai Makiyama, 1935; and M. arenaria corbuloides Comfort, 1938 (Coan, Scott & Bernard, 2000; Zhang, Xu & Liu, 2012). One of the most contentious questions pertains to the validity of M. japonica, which has been reported in the western Pacific (Habe, 1955, 1977; Fujie, 1957, 1962; MacNeil, 1965; Okutani, 2000), with possible relict populations or adventitious individuals occurring in the Bering Sea (MacNeil, 1965; Carlton, 1979). Soft-shell clams in the western Pacific identified as M. japonica possess rough commarginal wrinkles, relatively shorter and more rounded posterior end and more impressed pallial line in comparison to M. arenaria (Jay, 1856; MacNeil, 1965). These differences have led some malacologists to accept M. japonica as a valid species (Habe, 1955, 1977; Fujie, 1957, 1962; MacNeil, 1965; Okutani, 2000), whereas others have reported that M. arenaria and M. japonica could not be separated by the statistical analysis of morphological characteristics, and therefore, M. japonica was a junior synonym of M. arenaria (Nagao & Inoue, 1941; Bernard, 1979, 1983; Coan et al., 2000; Huber, 2010; Zhang et al., 2012). Confounding these taxonomic uncertainties are conflicting descriptions of important morphological characters, such as the chondrophore in the left valve (MacNeil, 1965; Bernard, 1979). Although several studies have examined genetic relationships between North American and European populations of M. arenaria (Strasser & Barber, 2009; Layton, Martel & Hebert, 2014; Barco et al., 2016; Cross et al., 2016; Lasota et al., 2016), none have conducted a rigorous genetic analysis of M. arenaria s.l. from across its entire range (i.e. including the western Pacific). Recent molecular analyses have identified varying levels of cryptic diversity in several cosmopolitan and wide-ranging molluscs. For example, many temperate North American populations of the tellinid clam Macoma balthica Linnaeus, 1758 are now identified as Macoma petalum Valenciennes, 1821 (Väinölä, 2003; Nikula, Strelkov & Väinölä, 2007); the widespread African freshwater oyster Etheria elliptica Lamarck, 1807 consists of at least three distinct sibling species (Elderkin et al., 2016); and the aeolidiid nudibranchs Aeolidia papillosa Linnaeus, 1761 and Spurilla neapolitana Delle Chiaje, 1841 now comprise complexes of four and five sibling species, respectively (Carmona et al., 2014; Kienberger et al., 2016). Molecular analyses can be particularly powerful when combined with morphological and other biological data (e.g. spermatozoan ultramorphology; Carstensen et al., 2009) within a broader integrative taxonomic framework (Padial et al., 2010; Schlick-Steiner et al., 2010). Spermatozoan ultramorphology appears to be highly conserved and has been used to resolve taxonomic and phylogenetic questions at the species level in many bivalve genera, including Solen (Hodgson, Devilliers & Bernard, 1987), Donax (Hodgson, Bernard & Vanderhorst, 1990; Introini, Passos & Recco-Pimentel, 2013), Crassostrea (Gwo, Liou & Cheng, 1996), Bathymodiolus (Eckelbarger & Young, 1999), Nutricola (Geraghty, Russell & Dollahon, 2008) and Brachidontes (Torroglosa & Gimenez, 2015). Fine organization of soft-shell clam spermatozoa from the western Pacific was studied in M. japonica (Drozdov, Sharina & Tyurin, 2009) and M. arenaria oonogai (Kim, Chung & Park, 2011). In both studies, spermatozoa were classified as ect-aquasperm with a conical curved head consisting of an acrosomal complex and nucleus, a midpiece with four mitochondria and a flagellum with a 9 × 2 + 2 organization of microtubules. However, no previous studies have described the morphology of M. arenaria spermatozoa from the greater Atlantic or preformed range-wide comparisons. In this study, we report the first rigorous range-wide investigation of the soft-shell clam, using an integrative taxonomic approach based on mitochondrial cytochrome c oxidase subunit I (COI), mitochondrial 16S rRNA (16S) and nuclear 28S rRNA (28S) gene data, and quantitative sperm ultramorphology and qualitative shell morphological comparisons. We also report the first detailed description of spermatozoan ultramorphology of M. arenaria s.s. and provide updated distributional data. MATERIAL AND METHODS Specimen collection For DNA analyses, we collected soft-shell clams from Yellow Sea, Qingdao, China (n = 8); Barents Sea, Murmansk, Russia (n = 3); Vostok Bay, Sea of Japan, Russia (n = 1); and Chesapeake Bay, Maryland, USA (n = 10) (Fig. 1). We also collected a congener, Mya uzenensis Nomura et Zinbo, 1937 from Vostok Bay (n = 1). The authors follow Lutaenko & Noseworthy (2012) and accept M. uzenensis as a valid species and not a junior synonym of M. pseudoarenaria. Voucher specimens were deposited at the Marine Biological Museum, Chinese Academy of Sciences (MBM), Qingdao, China, and the Smithsonian Institution National Museum of Natural History (USNM), Washington, DC, USA (Table 1). For spermatozoan ultramorphological analyses, three males in each of the Yellow Sea, Barents Sea and Sea of Japan sampling sites were collected at pre-spawning and spawning stages of the reproductive cycle (for staging details, see Porter, 1974; Dzyuba & Maslennikova, 1987). Voucher specimens were deposited at the Zoological Museum, Far Eastern Federal University, Vladivostok, Russia (ZMFU). In addition to the samples mentioned above, we examined 135 lots of soft-shell clams (photographing a sub-set of samples), which included individuals previously identified as M. arenaria and M. japonica, from historical MBM and ZMFU collections to identify possible morphological differences between soft-shell clams collected in the western Pacific and northwestern Asia/northern Europe (e.g. Barents Sea, White Sea, North Sea, Baltic Sea; Fig. 2). Table 1. Sampling information of specimens collected for this study and GenBank accession numbers of their sequences Species Voucher Location GenBank numbers COI 16S 28S Mya arenaria MBM229022 Murmansk, Barents Sea, Russia (Arctic) KX534201 KX534176 KX534188 Mya arenaria MBM229023 Murmansk, Barents Sea, Russia (Arctic) KX534202 KX534177 KX534189 Mya arenaria MBM229024 Murmansk, Barents Sea, Russia (Arctic) KX534203 KX534178 KX534190 Mya arenaria USNM:IZ:1286707 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959460 KT959520 – Mya arenaria USNM:IZ:1286708 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959399 KT959490 – Mya arenaria USNM:IZ:1286709 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959418 KT959498 – Mya arenaria USNM:IZ:1286710 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959394 KT959487 – Mya arenaria USNM:IZ:1286892 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906055 – – Mya arenaria USNM:IZ:1286895 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906064 – – Mya arenaria USNM:IZ:1286896 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905943 – – Mya arenaria USNM:IZ:1286897 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905735 – – Mya arenaria USNM:IZ:1287437 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906072 – – Mya arenaria USNM:IZ:1287540 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906038 – – Mya japonica MBM229013 Qingdao, Yellow Sea, China (NW Pacific) KX534192 KX534167 KX534180 Mya japonica MBM229014 Qingdao, Yellow Sea, China (NW Pacific) KX534193 KX534168 – Mya japonica MBM229015 Qingdao, Yellow Sea, China (NW Pacific) KX534194 KX534169 KX534181 Mya japonica MBM229016 Qingdao, Yellow Sea, China (NW Pacific) KX534195 KX534170 KX534182 Mya japonica MBM229017 Qingdao, Yellow Sea, China (NW Pacific) KX534196 KX534171 KX534183 Mya japonica MBM229018 Qingdao, Yellow Sea, China (NW Pacific) KX534197 KX534172 KX534184 Mya japonica MBM229019 Qingdao, Yellow Sea, China (NW Pacific) KX534198 KX534173 KX534185 Mya japonica MBM229020 Qingdao, Yellow Sea, China (NW Pacific) KX534199 KX534174 KX534186 Mya japonica MBM229021 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534200 KX534175 KX534187 Mya uzenensis MBM229025 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534204 KX534179 KX534191 Species Voucher Location GenBank numbers COI 16S 28S Mya arenaria MBM229022 Murmansk, Barents Sea, Russia (Arctic) KX534201 KX534176 KX534188 Mya arenaria MBM229023 Murmansk, Barents Sea, Russia (Arctic) KX534202 KX534177 KX534189 Mya arenaria MBM229024 Murmansk, Barents Sea, Russia (Arctic) KX534203 KX534178 KX534190 Mya arenaria USNM:IZ:1286707 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959460 KT959520 – Mya arenaria USNM:IZ:1286708 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959399 KT959490 – Mya arenaria USNM:IZ:1286709 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959418 KT959498 – Mya arenaria USNM:IZ:1286710 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959394 KT959487 – Mya arenaria USNM:IZ:1286892 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906055 – – Mya arenaria USNM:IZ:1286895 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906064 – – Mya arenaria USNM:IZ:1286896 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905943 – – Mya arenaria USNM:IZ:1286897 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905735 – – Mya arenaria USNM:IZ:1287437 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906072 – – Mya arenaria USNM:IZ:1287540 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906038 – – Mya japonica MBM229013 Qingdao, Yellow Sea, China (NW Pacific) KX534192 KX534167 KX534180 Mya japonica MBM229014 Qingdao, Yellow Sea, China (NW Pacific) KX534193 KX534168 – Mya japonica MBM229015 Qingdao, Yellow Sea, China (NW Pacific) KX534194 KX534169 KX534181 Mya japonica MBM229016 Qingdao, Yellow Sea, China (NW Pacific) KX534195 KX534170 KX534182 Mya japonica MBM229017 Qingdao, Yellow Sea, China (NW Pacific) KX534196 KX534171 KX534183 Mya japonica MBM229018 Qingdao, Yellow Sea, China (NW Pacific) KX534197 KX534172 KX534184 Mya japonica MBM229019 Qingdao, Yellow Sea, China (NW Pacific) KX534198 KX534173 KX534185 Mya japonica MBM229020 Qingdao, Yellow Sea, China (NW Pacific) KX534199 KX534174 KX534186 Mya japonica MBM229021 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534200 KX534175 KX534187 Mya uzenensis MBM229025 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534204 KX534179 KX534191 View Large Table 1. Sampling information of specimens collected for this study and GenBank accession numbers of their sequences Species Voucher Location GenBank numbers COI 16S 28S Mya arenaria MBM229022 Murmansk, Barents Sea, Russia (Arctic) KX534201 KX534176 KX534188 Mya arenaria MBM229023 Murmansk, Barents Sea, Russia (Arctic) KX534202 KX534177 KX534189 Mya arenaria MBM229024 Murmansk, Barents Sea, Russia (Arctic) KX534203 KX534178 KX534190 Mya arenaria USNM:IZ:1286707 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959460 KT959520 – Mya arenaria USNM:IZ:1286708 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959399 KT959490 – Mya arenaria USNM:IZ:1286709 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959418 KT959498 – Mya arenaria USNM:IZ:1286710 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959394 KT959487 – Mya arenaria USNM:IZ:1286892 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906055 – – Mya arenaria USNM:IZ:1286895 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906064 – – Mya arenaria USNM:IZ:1286896 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905943 – – Mya arenaria USNM:IZ:1286897 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905735 – – Mya arenaria USNM:IZ:1287437 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906072 – – Mya arenaria USNM:IZ:1287540 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906038 – – Mya japonica MBM229013 Qingdao, Yellow Sea, China (NW Pacific) KX534192 KX534167 KX534180 Mya japonica MBM229014 Qingdao, Yellow Sea, China (NW Pacific) KX534193 KX534168 – Mya japonica MBM229015 Qingdao, Yellow Sea, China (NW Pacific) KX534194 KX534169 KX534181 Mya japonica MBM229016 Qingdao, Yellow Sea, China (NW Pacific) KX534195 KX534170 KX534182 Mya japonica MBM229017 Qingdao, Yellow Sea, China (NW Pacific) KX534196 KX534171 KX534183 Mya japonica MBM229018 Qingdao, Yellow Sea, China (NW Pacific) KX534197 KX534172 KX534184 Mya japonica MBM229019 Qingdao, Yellow Sea, China (NW Pacific) KX534198 KX534173 KX534185 Mya japonica MBM229020 Qingdao, Yellow Sea, China (NW Pacific) KX534199 KX534174 KX534186 Mya japonica MBM229021 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534200 KX534175 KX534187 Mya uzenensis MBM229025 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534204 KX534179 KX534191 Species Voucher Location GenBank numbers COI 16S 28S Mya arenaria MBM229022 Murmansk, Barents Sea, Russia (Arctic) KX534201 KX534176 KX534188 Mya arenaria MBM229023 Murmansk, Barents Sea, Russia (Arctic) KX534202 KX534177 KX534189 Mya arenaria MBM229024 Murmansk, Barents Sea, Russia (Arctic) KX534203 KX534178 KX534190 Mya arenaria USNM:IZ:1286707 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959460 KT959520 – Mya arenaria USNM:IZ:1286708 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959399 KT959490 – Mya arenaria USNM:IZ:1286709 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959418 KT959498 – Mya arenaria USNM:IZ:1286710 Chesapeake Bay, Maryland, USA (NW Atlantic) KT959394 KT959487 – Mya arenaria USNM:IZ:1286892 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906055 – – Mya arenaria USNM:IZ:1286895 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906064 – – Mya arenaria USNM:IZ:1286896 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905943 – – Mya arenaria USNM:IZ:1286897 Chesapeake Bay, Maryland, USA (NW Atlantic) KU905735 – – Mya arenaria USNM:IZ:1287437 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906072 – – Mya arenaria USNM:IZ:1287540 Chesapeake Bay, Maryland, USA (NW Atlantic) KU906038 – – Mya japonica MBM229013 Qingdao, Yellow Sea, China (NW Pacific) KX534192 KX534167 KX534180 Mya japonica MBM229014 Qingdao, Yellow Sea, China (NW Pacific) KX534193 KX534168 – Mya japonica MBM229015 Qingdao, Yellow Sea, China (NW Pacific) KX534194 KX534169 KX534181 Mya japonica MBM229016 Qingdao, Yellow Sea, China (NW Pacific) KX534195 KX534170 KX534182 Mya japonica MBM229017 Qingdao, Yellow Sea, China (NW Pacific) KX534196 KX534171 KX534183 Mya japonica MBM229018 Qingdao, Yellow Sea, China (NW Pacific) KX534197 KX534172 KX534184 Mya japonica MBM229019 Qingdao, Yellow Sea, China (NW Pacific) KX534198 KX534173 KX534185 Mya japonica MBM229020 Qingdao, Yellow Sea, China (NW Pacific) KX534199 KX534174 KX534186 Mya japonica MBM229021 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534200 KX534175 KX534187 Mya uzenensis MBM229025 Vostok Bay, Peter the Great Bay, Sea of Japan, Russia (NW Pacific) KX534204 KX534179 KX534191 View Large Figure 1. View largeDownload slide Locations of Mya japonica (red triangles), Mya arenaria (red dots) and Mya uzenensis (red diamond) collected in this study, and other known localities of Mya japonica (yellow triangles), Mya arenaria (yellow dots) and Mya truncata (yellow squares) of which sequences were used in genetic analyses. Figure 1. View largeDownload slide Locations of Mya japonica (red triangles), Mya arenaria (red dots) and Mya uzenensis (red diamond) collected in this study, and other known localities of Mya japonica (yellow triangles), Mya arenaria (yellow dots) and Mya truncata (yellow squares) of which sequences were used in genetic analyses. Figure 2. View largeDownload slide Mya shell morphology. A–E, Mya japonica Jay, 1857, ZMFU no. 10069/Bv-514, sampled from the open part of Peter the Great Bay, Bolshoy Pelis Island, Russia. F–J, Mya japonica Jay, 1857, ZMFU no. 23301/Bv-3785, sampled from semi-enclosed part of Peter the Great Bay, Vityaz Bay, Russia. K–O, Mya japonica Jay, 1857, MBM229023, sampled from Qingdao, Yellow Sea, China. P–T, Mya arenaria Linnaeus, 1758, MBM229013, sampled from Murmansk, Barents Sea, Russia. Figure 2. View largeDownload slide Mya shell morphology. A–E, Mya japonica Jay, 1857, ZMFU no. 10069/Bv-514, sampled from the open part of Peter the Great Bay, Bolshoy Pelis Island, Russia. F–J, Mya japonica Jay, 1857, ZMFU no. 23301/Bv-3785, sampled from semi-enclosed part of Peter the Great Bay, Vityaz Bay, Russia. K–O, Mya japonica Jay, 1857, MBM229023, sampled from Qingdao, Yellow Sea, China. P–T, Mya arenaria Linnaeus, 1758, MBM229013, sampled from Murmansk, Barents Sea, Russia. Prior to analyses, all specimens (field collected and museum vouchers) were identified based on shell morphological differences of these two putative species (M. arenaria or M. japonica) recognized by early malacologists (Habe, 1955, 1977; Fujie, 1957, 1962; MacNeil, 1965). Specifically, we examined shell shape, noting profile of the posterior end, shell texture (regularity and depth of commarginal growth lines), the pallial line, and the shape, size, and orientation of the chondrophore in the left valve. Soft-shell clams possessing an elongated shell shape with even commarginal growth lines and lightly impressed pallial line were classified as M. arenaria, whereas soft-shell clams with a more truncate shell shape with rough commarginal wrinkles and more deeply impressed pallial line were classified as M. japonica. Molecular work DNA extraction, amplification and sequencing Thirteen specimens of M. arenaria, nine specimens of M. japonica and one specimen of M. uzenensis were collected for molecular analyses. We extracted total DNA from the samples using the Marine Animal Genomic DNA Extraction Kit (Tiangen Biotech, Beijing, China) according to manufacturer’s instructions. North American samples were sequenced following Aguilar et al. (2016) using the Primer sets jgLCO1490/jgHCO2198 (COI) and 16Sar-L/16Sb (16S). For Eurasian samples, fragments of other specimens of the mitochondrial COI gene were amplified using primers of Lco1490: 5′-GGT CAA CAA ATC ATA AAG ATA TTG G-3′, Hco2198: 5′-TAA ACT TCA GGG TGA CCA AAA AAT CA-3′ (Folmer et al., 1994). Fragments of the mitochondrial 16S rRNA gene were amplified using 16Sa: 5′-CGC CTG TTT ATC AAA AAC AT-3′ (Xiong & Kocher, 1991), 16Sb: 5′-CTC CGG TTT GAA CTC AGA TCA-3′ (Edgecombe, Giribet & Wheeler, 2002). The nuclear 28S rRNA gene was amplified using 28Sa: 5′-GAC CCG TCT TGA AAC ACG GA-3′, 28Sb: 5′-TCG GAA GGA ACC AGC TAC-3′ (Whiting et al., 1997). The 25-µL reaction for the amplification contained 12 µL 2×Es TaqMasterMix (CWBio Co., Ltd, Beijing, China), 2 µL of template DNA (50 ng/µL), 1 µL of each primer (10 M) and 9 µL dH2O. Amplification reactions were conducted in a T100 Thermal Cycler (Bio-Rad Laboratories, Inc.) with the following thermal profile: 95 °C for 3 min; 35 cycles of 95 °C for 30 s, 50–55 °C for 45 s, 72 °C for 1 min, and a further 10-min elongation at 72 °C. The PCR products were purified and sequenced by BGI Tech Solutions Co., Ltd. Species delimitation analyses To discriminate species, we used an integrative taxonomic approach including pairwise uncorrected p-distances, phylogenetic trees, Automatic Barcode Gap Discovery (ABGD) and spermatozoan ultramorphological and qualitative shell morphological data. In addition to sequences generated from our field collections, we included all Mya COI (> 500 bp), 16S and 28S gene sequences available in GenBank and BOLD (Ratnasingham & Hebert, 2007) in the generation of p-distances and in the ABGD analyses, most of which were associated with published articles (Supporting Information, Table S1). In total, 175 COI sequences (including four Potamocorbula sequences, which were originally submitted under the name Corbula amurensis in GenBank as an outgroup) were used for ABGD and p-distance analyses, and 26 16S sequences (including four outgroup Potamocorbula sequences) were also used for p-distance analyses. To improve clarity, sequences without associated publications or location information were omitted (n = 49) from phylogenetic trees. The ABGD (Puillandre et al., 2012) settings were as follows: steps = 10, X = 1, Pmin = 0.001, Pmax = 0.1, Nb bins = 20, and Jukes-Cantor (JC69) and Kimura (K80) TS/TV. The p-distances within and among each species grouping were estimated with MEGA v.6 (Tamura et al., 2013). Phylogenetic analyses Sequences of each region were aligned independently using MAFFT v.7 (Katoh & Standley, 2013) with the G-INS-i and Q-INS-i algorithms for the protein-coding and ribosomal regions. The alignments were further checked by eye. A combined analysis was conducted with concatenated sequences from the three genes. Three genes from the same individual were concatenated together as a sequence in SequenceMatrix v 1.8 (Vaidya, Lohman & Meier, 2011). Sequences from Potamocorbula amurensis (Schrenck, 1861) in the same superfamily, Myoidea, were selected as an outgroup. The best-fitted evolutionary models were selected by Akaike information criterion as implemented in jModeltest2 (Darriba et al., 2012). Maximum likelihood (ML) analysis was carried out using PhyML-3.1 (Guindon et al., 2010). Node support came from a majority-rule consensus tree of 1000 bootstrap replicates. For the ML bootstraps, we consider values < 70 as low, 70‒95 as moderate and ≥ 95 as high, following Hillis & Bull (1993). Bayesian inference (BI) analysis was carried out using MrBayes v3.2.6 (Ronquist et al., 2012) in CIPRES Science Gateway. We estimated posterior probability (PP) using four chains running 10000000 generations sampling every 100 generations. The first 25% of sampled trees were considered burn-in trees, and Tracer v1.6 was used to confirm that post-burn-in trees yielded an effective sample size of > 200 for all parameters. All remaining trees were used to calculate PPs using a 50% majority-rule consensus. For the Bayesian PPs, we considered values < 94% as low and ≥ 95% as high following Alfaro, Zoller & Lutzoni (2003). Divergence time estimates We used a Bayesian Markov chain Monte Carlo method implemented in Beast version 2.4.0 (Bouckaert et al., 2014) to estimate the evolutionary ages of internal nodes in the tree topology derived from the combined phylogenetic analyses. The time to most recent common ancestor (T), that is the time of divergence, was estimated as T = K/2r, where K is the K2P genetic distance, r is the nucleotide substitution rate. Substitution rates (r) were estimated as percentage per lineage per million years (Myr), so they equal one half the pairwise sequence divergence per Myr, that is divergence rate. Our estimates of T are based on the mitochondrial COI divergence between cognate species of Arcidae isolated across the Isthmus of Panama at a rate of 0.7–1.2%/Myr (Marko, 2002); thus, the substitution rate (r) equals 0.35–0.6%/Myr. A relaxed, uncorrelated lognormal clock with a general time reversible (GTR) substitution model was inferred for each partition, and a Yule speciation process was assumed for the tree prior. The analysis was run for 100 million generations, sampled every 1000 generations and a burn-in of 25%. Convergence and mixing were assessed using Tracer v. 1.6 by examining log-likelihood values across generations and ensuring that all remaining trees yielded an effective sample size for all parameters. Results were visualized using Tracer v. 1.6 and FigTree v. 1.4.2. Spermatozoan ultramorphology Transmission electron microscopes preparation Small pieces of gonadal tissues were pre-fixed in 2.5% (v/v) glutaraldehyde in 0.1 M cacodylate buffer (pH 7.5), with the addition of NaCl to provide appropriate tonicity of sea water, for 2 h at 4 °С. The specimens were then rinsed several times in the same buffer and held at 4 °C until analysis at the National Scientific Center of Marine Biology (Vladivostok, Russia). In the laboratory, gonadal tissue samples were post-fixed in 2% (w/v) osmium tetroxide in 0.1 M cacodylate buffer (pH 7.5) for 2 h in the dark at room temperature. The specimens were washed in distilled water and then dehydrated using a graded ethanol series and acetone. The dehydrated specimens were embedded in Araldite-EMbed-812 resin. Ultrathin sections (c. 50 nm) were obtained using a Leica UC6 ultramicrotome equipped with a diamond knife. The sections were stained with 2% aqueous uranyl acetate and Reynolds’ lead citrate (Reynolds, 1963) and were observed under Zeiss Libra 120 and Zeiss Libra 200FE transmission electron microscopes (Far East Center of Electron Microscopy, National Scientific Center of Marine Biology FEB RAS). Measurements of spermatozoa To investigate possible differences in spermatozoan ultramorphology among soft-shell clams, we examined three males (five spermatozoa per male) from each species/geographical groupings: M. arenaria from the Barents Sea, M. japonica from the Sea of Japan and M. japonica from the Yellow Sea. Only those spermatozoa in which the following compartments could be simultaneously observed in a longitudinal projection were measured: (1) acrosomal complex with cylindrical subacrosomal space, from the base of anterior nuclear fossa to the apex of acrosomal vesicle; (2) cone-shaped nucleus with the gradual dilatation from the apical to the basal part without any narrowing; and (3) midpiece with two centrioles. Because of the curved shape of the nucleus, we measured the maximal Feret’s diameter to estimate the nuclear length. We measured the length of the acrosomal complex from the invagination of the anterior nuclear fossa to the terminal point at the apical part of the acrosomal vesicle, and the width of the acrosomal complex at its base. All measurements were calculated with ImageJ version 1.49b (Schneider, Rasband & Eliceiri, 2012). To compare the size and shape of the acrosomal complex, we calculated volumes and length to width ratios. Based on the close-to-conical profile of the acrosomal complexes in the sections, we used the following formula to estimate volume: V = 1/3πR2H, where R is ½ of the width of the acrosomal complex and H is the length of the acrosomal complex. Statistical analyses We used one-way analysis of variance (ANOVA) to test for significant differences in mean nucleus length, mean acrosomal volume and mean acrosomal length to width ratio among the three species/geographical groupings: M. arenaria from the Barents Sea, M. japonica from the Sea of Japan and M. japonica from the Yellow Sea. Effect size (η2) for ANOVA was calculated using the corresponding function of R add-on package ‘lsr’ version 0.5 (Navarro, 2014). Multiple comparisons among geographical areas were made using Tukey’s honestly significant difference (Tukey’s HSD) tests. Prior to analysis, all data sets were tested for normality with Shapiro-Wilk tests and for homogeneity of variance using Levene’s tests and were determined to meet these assumptions of parametric analyses. Differences were considered statistically significant at P < 0.05 for all tests. We performed all statistical analyses using R version 3.2.4 (R Core Team, 2015) running on a Debian GNU/Linux workstation. RESULTS Morphology Shell characteristics and morphological identifications Based on our visual examinations, the samples from the western Pacific morphologically corresponded well with M. japonica and the northeast Pacific/Atlantic/Arctic soft-shell calms were morphologically identical with M. arenaria (Fig. 2). Detailed statistical analyses of morphological characteristics of soft-shell clams classified as M. arenaria and M. japonica have been reported previously (see Bernard, 1979), but here we describe minor differences in shell characteristics between them. In general, M. arenaria possessed a more elongated shell, with a longer and thinner posterior end and regular thin commarginal growth lines (Fig. 2P–T). Mya japonica was higher, with rough commarginal wrinkles and a relatively shorter and more rounded posterior end (Fig. 2A–O). In addition, the pallial line of M. arenaria was less impressed than in M. japonica. However, there was considerable variation in the shape, relative size and orientation of the left valve chondrophore among all the clams examined and not useful for discriminating between species. Two distinct morphological forms of M. japonica were observed in Peter the Great Bay. A thin, ‘typical’ form of M. japonica, which was elongated with regular growth lines, occurring mostly in semi-enclosed bays and other protected seashore areas, often with decreased salinity and a low energy environment (Fig. 2F–J). Conversely, a ‘thick-shelled’ form with irregular growth lines occurred in open, wave-influenced bays (Fig. 2A–E). Intermediate specimens were also found. Comparison of spermatozoan ultramorphology All spermatozoa examined in the present study were ect-aquasperm, a presumed primitive sperm type often associated with aquatic invertebrates with external fertilizing, free-swimming gametes (Eckelbarger & Young, 1999). Spermatozoa of all the specimens examined consisted of an acrosomal complex, nucleus, midpiece and flagellum. The data on dimensions of nuclei and acrosomal complexes are summarized in Table 2. The sperm were all rotationally asymmetrical due to the curvature of the nucleus or the position of the acrosomal complex (Fig. 3A–C). The midpieces of all spermatozoa examined were similar, containing four spherical mitochondria (Fig. 3G–I), and orthogonally arranged proximal and distal centrioles (Fig. 3J–L). In addition, the flagellum of all spermatozoa possessed a typically organized axoneme (9 × 2 + 2), posterior to the distal centriole. Table 2. Descriptive statistics of spermatozoan compartment dimensions in Mya arenaria sampled from the Barents Sea and in Mya japonica sampled from the Sea of Japan and the Yellow Sea (mean ± SD) Mya arenaria from the Barents Sea Mya japonica from the Sea of Japan Mya japonica from the Yellow Sea Nucleus length (µm) 1.98 ± 0.16 2.97 ± 0.13 3.08 ± 0.18 Acrosomal complex length (µm) 0.74 ± 0.06 0.80 ± 0.05 0.82 ± 0.04 Acrosomal complex width (µm) 0.61 ± 0.05 0.55 ± 0.04 0.57 ± 0.03 Acrosomal complex length to width ratio 1.22 ± 0.11 1.46 ± 0.08 1.46 ± 0.09 Acrosomal complex volume (µm3) 0.073 ± 0.015 0.064 ± 0.012 0.069 ± 0.010 Mya arenaria from the Barents Sea Mya japonica from the Sea of Japan Mya japonica from the Yellow Sea Nucleus length (µm) 1.98 ± 0.16 2.97 ± 0.13 3.08 ± 0.18 Acrosomal complex length (µm) 0.74 ± 0.06 0.80 ± 0.05 0.82 ± 0.04 Acrosomal complex width (µm) 0.61 ± 0.05 0.55 ± 0.04 0.57 ± 0.03 Acrosomal complex length to width ratio 1.22 ± 0.11 1.46 ± 0.08 1.46 ± 0.09 Acrosomal complex volume (µm3) 0.073 ± 0.015 0.064 ± 0.012 0.069 ± 0.010 Number of samples: 15 spermatozoa in each sampling site. View Large Table 2. Descriptive statistics of spermatozoan compartment dimensions in Mya arenaria sampled from the Barents Sea and in Mya japonica sampled from the Sea of Japan and the Yellow Sea (mean ± SD) Mya arenaria from the Barents Sea Mya japonica from the Sea of Japan Mya japonica from the Yellow Sea Nucleus length (µm) 1.98 ± 0.16 2.97 ± 0.13 3.08 ± 0.18 Acrosomal complex length (µm) 0.74 ± 0.06 0.80 ± 0.05 0.82 ± 0.04 Acrosomal complex width (µm) 0.61 ± 0.05 0.55 ± 0.04 0.57 ± 0.03 Acrosomal complex length to width ratio 1.22 ± 0.11 1.46 ± 0.08 1.46 ± 0.09 Acrosomal complex volume (µm3) 0.073 ± 0.015 0.064 ± 0.012 0.069 ± 0.010 Mya arenaria from the Barents Sea Mya japonica from the Sea of Japan Mya japonica from the Yellow Sea Nucleus length (µm) 1.98 ± 0.16 2.97 ± 0.13 3.08 ± 0.18 Acrosomal complex length (µm) 0.74 ± 0.06 0.80 ± 0.05 0.82 ± 0.04 Acrosomal complex width (µm) 0.61 ± 0.05 0.55 ± 0.04 0.57 ± 0.03 Acrosomal complex length to width ratio 1.22 ± 0.11 1.46 ± 0.08 1.46 ± 0.09 Acrosomal complex volume (µm3) 0.073 ± 0.015 0.064 ± 0.012 0.069 ± 0.010 Number of samples: 15 spermatozoa in each sampling site. View Large Figure 3. View largeDownload slide Sperm fine morphology in the studied Mya species/populations. A, D, G, J, Mya arenaria sampled from the Barents Sea. B, E, H, K, Mya japonica sampled from the Sea of Japan. C, F, I, L, Mya japonica sampled from the Yellow Sea. A–C, longitudinal section of spermatozoon consisted of acrosomal complex (ac), nucleus (n) with anterior (anf) and posterior nuclear fossae (pnf), midpiece (mp) and flagellum (f), scale bar = 1 µm. D–F, acrosomal complex included acrosomal vesicle (av) with an inner moderate electron-dense layer (asterisk) and outer electron-dense layer (arrow) and subacrosomal space (sas) filled with subacrosomal material, scale bar = 0.3 µm. G–I, transversal section of midpiece showing four mitochondria (m) around distal centriole (dc), scale bar = 0.5 µm. J–L, longitudinal section of midpiece. Proximal (pc) and distal centrioles (dc) located orthogonally to each other, scale bar = 0.5 µm. Figure 3. View largeDownload slide Sperm fine morphology in the studied Mya species/populations. A, D, G, J, Mya arenaria sampled from the Barents Sea. B, E, H, K, Mya japonica sampled from the Sea of Japan. C, F, I, L, Mya japonica sampled from the Yellow Sea. A–C, longitudinal section of spermatozoon consisted of acrosomal complex (ac), nucleus (n) with anterior (anf) and posterior nuclear fossae (pnf), midpiece (mp) and flagellum (f), scale bar = 1 µm. D–F, acrosomal complex included acrosomal vesicle (av) with an inner moderate electron-dense layer (asterisk) and outer electron-dense layer (arrow) and subacrosomal space (sas) filled with subacrosomal material, scale bar = 0.3 µm. G–I, transversal section of midpiece showing four mitochondria (m) around distal centriole (dc), scale bar = 0.5 µm. J–L, longitudinal section of midpiece. Proximal (pc) and distal centrioles (dc) located orthogonally to each other, scale bar = 0.5 µm. Nuclei were curved in all three geographical groups, but in M. arenaria from the Barents Sea, the curvature of the nuclei was less obvious than that in M. japonica from the Sea of Japan and the Yellow Sea (Fig. 3A–C). All nuclei possessed electron-dense chromatin with irregularly shaped electron-lucent lacunae. Anterior portions of the nuclei were slightly invaginated and formed an anterior nuclear fossa that was located subterminally (in M. arenaria) or apically (in M. japonica; Fig. 3A–C). Posterior portions had a well-defined invagination (posterior nuclear fossa) where a proximal centriole was located, and indentations were in the regions of contact of the nucleus and mitochondria (Fig. 3A–C, J–L). There were statistically significant differences in mean nuclei length (ANOVA: F(2,42) = 219.51, P < 0.05, η2 = 0.91). The Tukey’s HSD test indicated that M. arenaria (Barents Sea) had significantly shorter nuclei compared to M. japonica from the Sea of Japan and Yellow Sea (P < 0.05; Fig. 4A). No statistically significant difference in nuclei length between M. japonica collected from the Sea of Japan and Yellow Sea was found (P = 0.165). Figure 4. View largeDownload slide Whisker plots showing results of Tukey’s honestly significant difference test for difference in mean nucleus lengths (A) and acrosomal L/W ratios (B) among groups of Mya. Black dots are the difference in pairwise comparisons between sampling locations of the means and the whiskers are 95% confidence intervals of these differences. Figure 4. View largeDownload slide Whisker plots showing results of Tukey’s honestly significant difference test for difference in mean nucleus lengths (A) and acrosomal L/W ratios (B) among groups of Mya. Black dots are the difference in pairwise comparisons between sampling locations of the means and the whiskers are 95% confidence intervals of these differences. In all spermatozoa examined, the acrosomal complexes were cone-shaped and consisted of a membrane-bounded acrosomal vesicle and subacrosomal space (Fig. 3D–F). Within the acrosomal vesicle, two regions with different electron density were distinguished: the inner with moderate electron density and the outer with high electron density (Fig. 3D–F). The acrosomal vesicle had a deep posterior invagination that extended to its apical portion and formed a subacrosomal space filled with flake-like subacrosomal material of low electron density. The distance between the acrosomal complex and the nucleus was narrow in M. arenaria (Barents Sea), but wider M. japonica (Sea of Japan and Yellow Sea; Fig. 3D–F). This is most likely the reason why the acrosomal complex lengths were generally larger relative to width in M. japonica (Sea of Japan and Yellow Sea) than M. arenaria (Barents Sea; Table 2). There was a statistically significant difference in mean acrosomal length to width ratios (ANOVA: F(2,42) = 30.83, P < 0.05, η2 = 0.60). A Tukey’s HSD test indicated that acrosomal complexes of the spermatozoa in M. arenaria (Barents Sea) were shorter relative to width than those in M. japonica (Sea of Japan and Yellow Sea; Table 2; P < 0.05), while no statistically significant difference in length to width ratios between Mya from the Sea of Japan and Yellow Sea was found (P = 1; Fig. 4B). There were no significant differences in mean acrosomal complex volumes among the different geographical regions (ANOVA: F(2,42) = 1.94, P = 0.156, η2 = 0.08). Molecular analyses Barcoding analyses In total, 52 new sequences were obtained and deposited in GenBank (Table 1). The alignment data set comprised 637 nucleotide positions for COI, 496 positions for 16S and 297 positions for 28S. The concatenated data set yielded a sequence alignment of 1430 positions. The ABGD analysis based on the COI alignment produced five distinct groupings, with P values ranging from 0.022 to 0.100 based on both the Jukes-Cantor (JC69) and Kimura (K80) TS/TV models (Table 3). These groupings corresponded with the five clades/subclades in the COI trees (Table 3; Fig. 5B–D). The COI p-distances within the former four groups were low (0–2.39%), while the p-distances among the five groups ranged from 7.17 to 19.12% (Table 3; Supporting Information, Table S2). The minimum interspecific p-distance among these groups was much higher than the maximum one within ABGD Group 1 (7.17 vs. 2.39%; Table 3), which indicated all the western Pacific soft-shell clams (those morphologically identified M. japonica by the authors and those classified as M. arenaria – the valid species name at the time of identification), plus two individuals from British Columbia (DFO097-11 + RBCMI219-14) represent a single species (i.e. M. japonica). After reassigning all western Pacific M. arenaria and DFO097-11 + RBCMI219-14 as M. japonica, the COI intraspecific p-distances within Mya were 0–7.57%, while the interspecific p-distances among Mya species ranged from 11.16 to 19.12% (Fig. 5A; Table 3; Supporting Information, Table S2), indicating a barcoding gap among Mya species (Fig. 5A). However, if the two Mya truncata clades (ABGD Groups 3 and 4) are considered distinct species, intraspecific p-distances decrease to 0–2.39%, while interspecific distances range from 7.17 to 19.12%. For 16S, the intraspecific p-distances ranged from 0 to 0.97% within M. arenaria and M. japonica, while the interspecific p-distances ranged from 5.3 to 13.80% (Supporting Information, Table S3), which also indicates a barcoding gap between M. arenaria and M. japonica. Table 3. Uncorrected pairwise distances at COI within and between the ABGD groups Group 1 Group 2 Group 3 Group 4 Group 5 Group 1 All Mya japonica + Mya arenaria QWEAS092- 15, QWEAS093-15, DFO097-11, RBCMI219- 14, KP976290, KP976289, KP976288, KP976287, KP976286, KJ125420, KJ125421, JQ267795 0‒2.39% Group 2 Remaining Mya arenaria (142 sequences) 11.55–13.94% 0‒1.99% Group 3 Mya truncata KF644116, KF644129 17.53–18.33% 15.54–17.53% 0.40% Group 4 Mya truncata KF644154, KF643769, KF643675, KF643403 17.53–18.73% 17.53–19.12% 7.17–11.55% 0 Group 5 Mya uzenensis KX534204 17.53–18.33% 16.33–17.93% 11.16–11.55% 13.15% ‒ Group 1 Group 2 Group 3 Group 4 Group 5 Group 1 All Mya japonica + Mya arenaria QWEAS092- 15, QWEAS093-15, DFO097-11, RBCMI219- 14, KP976290, KP976289, KP976288, KP976287, KP976286, KJ125420, KJ125421, JQ267795 0‒2.39% Group 2 Remaining Mya arenaria (142 sequences) 11.55–13.94% 0‒1.99% Group 3 Mya truncata KF644116, KF644129 17.53–18.33% 15.54–17.53% 0.40% Group 4 Mya truncata KF644154, KF643769, KF643675, KF643403 17.53–18.73% 17.53–19.12% 7.17–11.55% 0 Group 5 Mya uzenensis KX534204 17.53–18.33% 16.33–17.93% 11.16–11.55% 13.15% ‒ View Large Table 3. Uncorrected pairwise distances at COI within and between the ABGD groups Group 1 Group 2 Group 3 Group 4 Group 5 Group 1 All Mya japonica + Mya arenaria QWEAS092- 15, QWEAS093-15, DFO097-11, RBCMI219- 14, KP976290, KP976289, KP976288, KP976287, KP976286, KJ125420, KJ125421, JQ267795 0‒2.39% Group 2 Remaining Mya arenaria (142 sequences) 11.55–13.94% 0‒1.99% Group 3 Mya truncata KF644116, KF644129 17.53–18.33% 15.54–17.53% 0.40% Group 4 Mya truncata KF644154, KF643769, KF643675, KF643403 17.53–18.73% 17.53–19.12% 7.17–11.55% 0 Group 5 Mya uzenensis KX534204 17.53–18.33% 16.33–17.93% 11.16–11.55% 13.15% ‒ Group 1 Group 2 Group 3 Group 4 Group 5 Group 1 All Mya japonica + Mya arenaria QWEAS092- 15, QWEAS093-15, DFO097-11, RBCMI219- 14, KP976290, KP976289, KP976288, KP976287, KP976286, KJ125420, KJ125421, JQ267795 0‒2.39% Group 2 Remaining Mya arenaria (142 sequences) 11.55–13.94% 0‒1.99% Group 3 Mya truncata KF644116, KF644129 17.53–18.33% 15.54–17.53% 0.40% Group 4 Mya truncata KF644154, KF643769, KF643675, KF643403 17.53–18.73% 17.53–19.12% 7.17–11.55% 0 Group 5 Mya uzenensis KX534204 17.53–18.33% 16.33–17.93% 11.16–11.55% 13.15% ‒ View Large Figure 5. View largeDownload slide Genetic distance range (A), phylogenetic trees (B, D) and candidate species inferred from Automatic Barcode Gap Discovery (ABGD) (C). A, the COI intraspecific distance values with M. truncata are 0, 0.33, 6.93 and 7.26% marked with green spots. B, phylogenetic tree inferred by Bayesian analysis (BI). C, ABGD, based on the COI alignment, with both models of Jukes-Cantor (JC69) and Kimura (K80) TS/TV. D, phylogenetic tree inferred by maximum likelihood (ML). Figure 5. View largeDownload slide Genetic distance range (A), phylogenetic trees (B, D) and candidate species inferred from Automatic Barcode Gap Discovery (ABGD) (C). A, the COI intraspecific distance values with M. truncata are 0, 0.33, 6.93 and 7.26% marked with green spots. B, phylogenetic tree inferred by Bayesian analysis (BI). C, ABGD, based on the COI alignment, with both models of Jukes-Cantor (JC69) and Kimura (K80) TS/TV. D, phylogenetic tree inferred by maximum likelihood (ML). Phylogenetic analyses For COI, 16S rDNA, 28S rDNA and the concatenated alignments, the TIM3 + I + G, TPM1uf + I, HKY and TVM + G evolutionary models were the best-fitted models, respectively. We report only the results of trees generated from COI sequences (Fig. 5B, D), which better resolved relationships at species level, and trees from the concatenated regions (Fig. 6). Trees constructed from 28S were uninformative with low support values at all levels and were omitted. Trees constructed from 16S were nearly identical in topology with the concatenated regions (except for differences in support values) and were also omitted. The resulting two consensus trees generated using BI and ML analyses were generally consistent; thus, a single topology was presented for the concatenated regions with support values indicated on branches (Fig. 6). This topology was in agreement with the COI trees and showed species-level groupings with high support (Fig. 6). Figure 6. View largeDownload slide Bayesian (BI) and maximum likelihood (ML) consensus tree inferred from concatenated sequences of COI, 16S and 28S genes. Figure 6. View largeDownload slide Bayesian (BI) and maximum likelihood (ML) consensus tree inferred from concatenated sequences of COI, 16S and 28S genes. In all COI trees, the ABGD Group 1 (all western Pacific soft-shell clams and two individuals of M. arenaria from the northeast Pacific) formed a clade with high node support (BI 1.00, ML 100%) and the ABGD Group 2 (remaining M. arenaria) grouped together with low to high node support before clustering to other species (BI 1.00, ML 42%, Fig. 5B, D). In the COI BI tree and the concatenated trees, M. arenaria clustered with M. japonica with high support, while M. truncata formed a sister clade with M. uzenensis with full node support (Figs 5B, 6). In the COI BI tree, M. japonica and M. arenaria QWEAS092-15 + QWEAS093-15 + DFO097-11 + RBCMI219-14 branched early within Mya with high node support (ML 95%), while the remaining M. arenaria clustered with the clade M. truncata + M. uzenensis with low support (ML 25%, Fig. 5D). In all COI trees, M. truncata were separated into two subclades, KF644116 + KF644129 and KF644154 + KF643769 + KF643675 + KF643403 (BI 0.86, ML 88). Divergence time estimates Assuming r1 = 0.35%, M. japonica would have split from M. arenaria ~9.4 Myr [95% highest posterior density (HPD) bounds: 6.7–12.5 Myr]. Mya uzenensis and M. truncata separated 8.1 Myr (95% HPD bounds: 4.0–11.9 Myr; Fig. 7). Time (T) to the most recent common ancestor for the Mya would be 9.9–20.0 Myr. Alternatively, assuming a faster substitution rate (r2 = 0.6% per lineage per Myr), M. japonica would have split from its congener, M. arenaria, ~5.7 Myr (95% HPD bounds: 4.1–7.5 Myr). Mya uzenensis and M. truncata separated 4.7 Mya (95% HPD bounds: 2.4–6.8 Myr). T for Mya would be 5.8–11.4 Myr (Fig. 7). Figure 7. View largeDownload slide Phylogenetic relationships among bivalve molluscs in Mya with estimates of divergence times, based on COI sequences. Bars in the tree correspond to 95% credibility intervals. Values above the bars were estimated by substitution rate r1 = 0.35% and below by substitution rate r2 = 0.6%. Scar bars on the bottom indicate time before present (Myr) derived from the two substitution rates (r1 and r2). Figure 7. View largeDownload slide Phylogenetic relationships among bivalve molluscs in Mya with estimates of divergence times, based on COI sequences. Bars in the tree correspond to 95% credibility intervals. Values above the bars were estimated by substitution rate r1 = 0.35% and below by substitution rate r2 = 0.6%. Scar bars on the bottom indicate time before present (Myr) derived from the two substitution rates (r1 and r2). DISCUSSION Mya japonica is a valid species Although the morphological differences between the clams identified as M. arenaria and M. japonica in the present study were minor, the ABGD delimitation, reciprocal monophyly, genetic divergence and spermatozoan ultramorphology of soft-shell clams from the western Pacific strongly support the recognition of M. japonica as a valid species. Within Mya, we found two clades, one comprising M. arenaria and M. japonica and the other containing both other Mya species examined. While M. arenaria and M. japonica were most related to each other, there was deep divergence between M. arenaria from the Northeast Pacific, Atlantic, Mediterranean and Barents Sea (Arctic) and M. japonica from the Pacific, which was further supported by ABGD. Furthermore, we found significant differences in spermatozoan ultramorphology between M. arenaria (Barents Sea) and M. japonica (Sea of Japan and Yellow Sea). The spermatozoan differences and genetic divergence reported in this study suggest long-term reproductive isolation between these two morphologically similar species. The fine morphology of spermatozoa from M. japonica examined in this study was very similar to M. japonica from Peter the Great Bay, Sea of Japan (Drozdov et al., 2009) and ‘M. arenaria oonogai’ from the coastal Korea Strait of Samcheonpo, Korea (Kim et al., 2011), which should all represent M. japonica s.s. Similarly, previous junior synonyms of M. arenaria described from the western Pacific, such as M. oonogai and M. arenaria oonogai, probably represent M. japonica. Shell morphological examinations and intraspecific variation of M. japonica Although individuals identified as M. arenaria and M. japonica were similar in appearance, there were recognizable differences in general shell shape, shell texture and appearance of the pallial line. Although the shape, size and orientation of the left valve chondrophore are often important diagnostic characters among Mya and related bivalves (Bernard, 1979), we found considerable variation in the characteristics of the chondrophore among all the material examined and it was not useful in discriminating between species. Given the comparative nature of the characteristics reported here (e.g. less/more impressed, elongated/shorter etc.) and overall degree of intraspecific variation, it may be difficult to identify individual clams using the descriptions contained herein, particularly for those who have not had the benefit of studying many specimens. Thus, we strongly recommend molecular or other analyses to confidently separate M. arenaria and M. japonica. Makiyama (1934, 1935) reported two forms of soft-shell clam in Japan, a northern form described as M. japonica and the southern form described as M. oonogai. The latter was referred to as M. japonica oonogai by Habe (1955) and Fujie (1957, 1962) and M. arenaria oonogai by Kwon et al. (2001), Kwon, Park & Lee (1993) and Yoo (1967). Fujie (1957) not only agreed that northern forms were M. japonica but also described a forma ‘α’, which was shorter and more ovate. Nagao & Inoue (1941) described differences between northern and southern Japanese forms and noted many intermediate specimens. The northern group, which possessed a thick, relatively short and nearly equilateral shell that was somewhat well rounded or slightly truncated posteriorly, was considered M. arenaria (due to synonymy with M. japonica). The southern group, which was characterized by its generally thin, long and somewhat inequilateral shell with a narrowly rounded or frequently acuminated posterior extremity, was considered an unnamed species. Nagao & Inoue (1941) also reported that the chondrophore in the northern group was nearly perpendicular to the antero-posterior diameter of the shell, and was narrower and more deeply excavated, with its inner border less convex than in the southern group. The posterior ridge was also not well differentiated from the chondrophore. Similarly, Scarlato (1981) noted two morphological forms of M. japonica in Russian Far Eastern seas. One form is distributed in low-boreal (temperate) waters (Aniva Bay in southern Sakhalin, southern Kurile Islands and Peter the Great Bay), where this species reaches its maximum size and possesses the typical shape with a smooth surface and fine, regular growth lines. The other form was found in the northern part of its distributional range (Bering and Chukotsk seas, eastern Kamchatka and northern Sea of Okhotsk), where M. japonica never reaches its maximum length and often has an atypical shape, sometimes appearing misshapen or deformed. However, our examination of M. japonica from Peter the Great Bay indicates that these variations also occur in the same geographical locality and are related to different environments (Fig. 2) along with intermediate specimens. The co-occurrence of different morphological forms and intermediate specimens in a single geographical region strongly suggests that these forms are ecophenotypes. It should be noted that the ‘typical’ thin and regular form reported here fully complies with ‘M. japonica oonogai’ collected from Hokkaido and Nagasaki, Japan, and our ‘thick-shelled’ form is very similar to M. japonica forma ‘α’, as illustrated by Fujie (1957). In addition, the M. japonica of Fujie (1957) is intermediate between the two forms we reported here. The results of our phylogenetic and morphological analyses indicate that M. japonica is a single species, which displays a high degree of phenotypic plasticity in shell morphology, and many previously reported forms/species combinations probably reflect ecophenotypic variation or represent natural intraspecific variation. Distribution of M. arenaria and M. japonica The soft-shell clam M. arenaria is widely accepted to have a circumboreal distribution (Bernard, 1979, 1983; Strasser, 1999; Coan et al., 2000; Huber, 2010; Zhang et al., 2012). Previous molecular studies indicated a complex colonization history between Europe and North America, but these analyses lacked samples from the northwest Pacific – the putative range of M. japonica (Strasser & Barber, 2009; Layton et al., 2014; Barco et al., 2016; Cross et al., 2016; Lasota et al., 2016). Our genetic and spermatozoan ultramorphological analyses strongly suggest that soft-shell clams in the northwest Pacific represent M. japonica. Two publicly available sequences identified as ‘M. arenaria’ (the valid species name at the time of identification) from China (QWEAS092-15 + QWEAS093-15) were highly similar (< 2% difference COI; Supporting Information, Table S2) to western Pacific M. japonica. These individuals are probably M. japonica and do not reflect the presence of M. arenaria in the western Pacific. MacNeil (1965) reported that M. arenaria occurred in Japan only in the middle Miocene and that later and recent specimens from the Okhotsk Sea, west coast of Japan and China are M. japonica, which was supported by Lutaenko & Noseworthy (2012) and Scarlato (1981). Numerous other authors have also reported M. japonica from China, Japan, Korea and the Russian Far Eastern seas (Jay, 1856; Makiyama, 1935; Habe, 1955, 1977; Tchang, Qi & Li, 1955; Fujie, 1957; Scarlato, 1981; Okutani, 2000). While both M. arenaria and M. japonica were present along western North America at one time, there is no palaeontological or archaeological evidence that either species survived the Pleistocene glaciation, with the possible exception of relict populations in the Bering Sea (Carlton, 1979). Mya arenaria was first introduced into the northeastern Pacific via Eastern Oyster Crassostrea virginica Gmelin, 1791 seed plantings in San Francisco Bay, California, during the 1860–1870s (Carlton, 1979). Since that time, it has spread northward, by either intentional plantings or natural dispersal, to southern Alaska (Carlton, 1979; Strasser, 1999; Powers et al., 2006). In contrast, until the present study, there were no recent eastern Pacific records of M. japonica or any evidence of intentional or accidental introductions of western Pacific Mya into the eastern Pacific. Two publicly available sequences identified as ‘M. arenaria’ from the British Columbia, one collected from Quascilla Bay (DFO097-11) and the other from Graham Island (RBCMI219-14), were highly similar (< 2% difference COI; Supporting Information, Table S2) to all western Pacific M. japonica. These individuals, both collected in June 2011, are probably M. japonica and went unnoticed because of morphological similarities to M. arenaria, the expected species in British Columbia. Thus, these records represent the first modern documented occurrence of M. japonica in the northeast Pacific. In addition, after examining photographs of these two clams, the authors are confident in classifying them as M. japonica based on the morphological criteria used in this study. We are uncertain of the source of introduction or status of possible eastern Pacific M. japonica populations. However, similar to M. arenaria, M. japonica was probably introduced via oyster plantings, but in this case, the Pacific Oyster Crassostrea gigas Thunberg, 1793. These plantings, which occurred from the early 1910s until after the Second World War (Carlton, 1979; Lavoie, 2005), are also thought to have led to other bivalve introductions, such as the Japanese Littleneck Ruditapes philippinarum Adams & Reeve, 1850 and Quadrate Trapezium Neotrapezium liratum Reeve, 1843. Alternatively, M. japonica may have been introduced by ballast water, as is speculated with the Asian Semele Theora lubrica Gould, 1861 and Asian brackish-water clam P. amurensis in California (Fofonoff et al., 2003). In either case, M. japonica was probably introduced decades after M. arenaria became widespread and well established in the eastern Pacific (Carlton, 1979), which reduced the likelihood of detection of M. japonica in areas where M. arenaria was already present or in nearby localities. Furthermore, both individuals were larger-sized adults collected nearly 400 km apart, which suggests there may be viable populations of M. japonica in at least British Columbia or possibly other Pacific Northeast locations. Albeit unlikely, the presence of individuals with high COI similarly to M. japonica could represent F1 hybrids with maternal M. japonica contribution or possibly reflect historical or cryptic hybridization (Pfenninger, Reinhardt & Streit, 2002; Rees, Dioli & Kirkendall, 2003) between M. arenaria and M. japonica. While there were significant differences in sperm morphology between M. arenaria and M. japonica, we are uncertain if these differences (or ecological/behavioural differences) would prohibit hybridization. According to Laursen (1966), historical reports of M. arenaria from the Arctic Ocean are erroneous and should represent M. truncata ovata Jensen, 1900, which is currently recognized as a junior synonym of M. pseudoarenaria (see Huber, 2010). However, we identified M. arenaria from the Russian waters of the Barents Sea, a marginal sea of the Arctic Ocean, and it has also been reported further west near Forsøl, Norway (Crocetta & Turolla, 2011). A rigorous genetic analysis (incorporating both nuclear and mitochondrial markers) of M. arenaria and M. japonica collected from across the eastern and western Pacific, Bering Sea, Arctic Ocean and greater Atlantic is required to fully resolve the zoogeography of both species, determine the extent of M. japonica introductions in the eastern Pacific, identify any possible cryptic introductions of M. arenaria in the western Pacific, determine the identity of soft-shell clams in the Bering Sea and examine possible hybridization of M. arenaria and M. japonica. Divergence time and evolutionary/migration history We estimate that M. japonica diverged from M. arenaria 4.1–12.5 Myr; this agrees with the evolution of M. arenaria, which evolved during early Pliocene to late Miocene from its ancestor Mya fujiei MacNeil, 1965, which arose during middle Miocene (MacNeil, 1965; Strauch, 1972). Mya arenaria first appeared in the fossil record from Japan in late Miocene formations and almost at the same time from the eastern Pacific in California (Strauch, 1972). Mya japonica was also probably present in the eastern Pacific at least during the Pleistocene (MacNeil, 1965; Carlton, 1979). Although the formation of the Bering Strait during the Pliocene overlapped with the divergence of M. arenaria and M. japonica, it is likely that only M. arenaria crossed the Arctic to reach the western and eastern Atlantic (MacNeil, 1965; Bernard, 1979; Vermeij, 1989). Strauch (1972) suggested M. arenaria may also have migrated through the Central American Passage from the Pacific to the Atlantic, where it then spread to Europe in the late Pliocene (Strasser, 1999), but fossil evidence for this migratory route is lacking (Bernard, 1979). During the Pleistocene glaciation, M. arenaria was extirpated from all areas except the western Atlantic and M. japonica was extirpated from all areas except the western Pacific (Strasser, 1999). However, the status of fossil, archaeological and recent soft-shell clams in the Bering Sea requires further investigation, which could represent relict populations of M. arenaria or M. japonica, adventitious individuals of M. japonica from the western Pacific or an undescribed cryptic species (MacNeil, 1965; Bernard, 1979; Carlton, 1979). The evolutionary history of Mya is complex and there is much uncertainty regarding the identification and evolutionary relationships of fossil material (MacNeil, 1965) and their connection with extant forms, particularly as it relates to the synonymy of extant species and species identified solely from fossil material (Petersen, 1999). Due to the morphological similarities between M. arenaria and M. japonica and levels of phenotypic plasticity, it may be difficult to fully resolve their evolutionary and migration history. MacNeil (1965) and Strauch (1972) reported that the oldest records of M. truncata and Mya cuneiformis Böhm 1915 were from the middle Miocene, both of which were direct descendants of the late Oligocene or early Miocene species, Mya salmonensis Clark, 1932, while M. pseudoarenaria and Mya priapus Tilesius, 1822 derived from M. cuneiformis in the late Miocene. Nakashima (1999) determined that the oldest records of M. pseudoarenaria, M. cuneiformis and M. truncata were from the early Miocene and M. uzenensis may have directly evolved from M. salmonensis. Our results indicate that the separation time of M. uzenensis and M. truncata (3.9–21.0 Myr) is nearly the same as M. japonica and M. arenaria and that M. uzenensis and M. truncata share a common ancestor, but it is unclear if M. uzenensis directly evolved from M. salmonensis or M. cuneiformis. Mya truncata may reflect a species complex The ABGD analysis indicated that M. truncata was a candidate species complex (Fig. 5C), which included two distinct groups (ABGD Groups 3 and 4; Table 3). The COI genetic distances between individuals of these groups ranged from 6.93 to 7.26% (Table 3). While these are intermediate values between the maximum intraspecific p-distance (1.70%) within M. japonica + M. arenaria and the minimum interspecific p-distance (10.68%) among the other Mya species, this may reflect a true species-level difference, which is consistent with Petersen (1999) and Layton et al. (2014). At present, we are unable to fully resolve the taxonomic status of M. truncata. All the M. truncata COI sequences examined in this study were obtained from individuals collected in Nunavut, Canada (Layton et al., 2014). However, individuals in ABGD Group 3 were collected near Cornwallis Island, whereas individuals in ABGD Group 4 were collected around Igloolik Island, indicating possible spatial structure to the putative species complex. Specimens collected from the entire range of M. truncata (eastern Pacific, Arctic and northern Atlantic) and additional data (e.g. spermatozoan ultramorphology) are needed to clarify the taxonomic status of M. truncata. CONCLUSIONS Our results indicate that (1) M. japonica is a valid species, distributed in the northern Pacific, with introduced populations of unknown size and distribution in the northeast Pacific; (2) M. arenaria is distributed in the northeast Pacific, northern Atlantic, Mediterranean and Barents Sea (Arctic); (3) there were significant species-level differences in the spermatozoan ultramorphology of M. arenaria and M. japonica; (4) differences in shell morphology were minor, and given the levels of intraspecific variation noted among varying habitats, shell morphology alone may not be an effective method to differentiate M. arenaria and M. japonica; (5) the congeners, M. arenaria, M. japonica, M. uzenensis and M. truncata were all differentiated by DNA barcoding, with evidence of cryptic diversity in M. truncata from the Canadian Arctic; and (6) the timing of divergence between M. arenaria and M. japonica supports the evolutionary hypothesis that M. arenaria evolved during early Pliocene to late Miocene. Further collections of Mya species from across their entire range (especially Asia and Arctic Ocean/Bering Sea where collections have been scarce) and a rigorous analysis incorporating an integrative taxonomic approach using ecological, molecular, morphological and/or other biological data are required to fully resolve taxonomic relationships and distributions of these highly morpho-variable species. SUPPORTING INFORMATION Supporting information may be found in the online version of this article at the publisher’s web-site: Table S1. List of all specimens used for molecular analyses. Table S2. Interspecific and intraspecific uncorrected pairwise distances at COI among species of Mya and Potamocorbula. Table S3. Interspecific and intraspecific uncorrected pairwise distances at 16S among species of Mya and Potamocorbula. ACKNOWLEDGEMENTS This study was supported by the Special Funds for the Young Scholars of Taxonomy of the Chinese Academy of Sciences (No. ZSBR-009), the Bureau of International Cooperation Chinese Academy of Sciences, the Far Eastern Branch of the Russian Academy of Sciences (No. 15-I-6–007 o) and the Russian Foundation for Basic Research (No. RFFI-ofi 15-29-0245615). The authors are grateful to Amy Driskell (Smithsonian Laboratories of Analytical Biology) for sequencing the Chesapeake Bay Mya arenaria samples. We thank Prof. Jim Carlton (Williams College, USA) for helpful suggestions on the manuscript. Irina E. Volvenko (Zoological Museum, Far Eastern Federal University, Russia) kindly took photographs of shells. REFERENCES Aguilar R , Ogburn MB , Driskell AC , Weigt LA , Groves MC , Hines AH . 2016 . Gutsy genetics: identification of digested piscine prey items in the stomach contents of sympatric native and introduced warmwater catfishes via DNA barcoding . Environmental Biology of Fishes 100 : 325 – 336 . Google Scholar CrossRef Search ADS Alfaro ME , Zoller S , Lutzoni F . 2003 . Bayes or bootstrap? A simulation study comparing the performance of Bayesian Markov chain Monte Carlo sampling and bootstrapping in assessing phylogenetic confidence . Molecular Biology and Evolution 20 : 255 – 266 . Google Scholar CrossRef Search ADS PubMed Barco A , Raupach MJ , Laakmann S , Neumann H , Knebelsberger T . 2016 . Identification of North Sea molluscs with DNA barcoding . Molecular Ecology Resources 16 : 288 – 297 . Google Scholar CrossRef Search ADS PubMed Bernard FR . 1979 . Identification of the living Mya (Bivalvia: Myoida) . Venus: Japanese Journal of Malacology 38 : 185 – 204 . Bernard FR . 1983 . Catalogue of the living Bivalvia of the eastern Pacific Ocean: Bering Strait to Cape Horn . Ottawa : Love Printing Services Ltd . Bouchet P , Gofas S . 2013 . Mya Linnaeus, 1758 . In: MolluscaBase (2017) . World Register of Marine Species . Available at: http://marinespecies.org/aphia.php?p=taxdetails&id=138211 (accessed 10 May 2017 ). Bouckaert R , Heled J , Kühnert D , Vaughan T , Wu CH , Xie D , Suchard MA , Rambaut A , Drummond AJ . 2014 . BEAST 2: a software platform for Bayesian evolutionary analysis . PLoS Computational Biology 10 : e1003537 . Google Scholar CrossRef Search ADS PubMed Carlton JT . 1979 . History, biogeography, and ecology of the introduced marine and estuarine invertebrates of the Pacific coast of North America . Unpublished Ph.D. Thesis, University of California . Carmona L , Lei BR , Pola M , Gosliner TM , Valdés Á , Cervera JL . 2014 . Untangling the Spurilla neapolitana (Delle Chiaje, 1841) species complex: a review of the genus Spurilla Bergh, 1864 (Mollusca: Nudibranchia: Aeolidiidae) . Zoological Journal of the Linnean Society 170 : 132 – 154 . Google Scholar CrossRef Search ADS Carstensen D , Laudien J , Leese F , Arntz W , Held C . 2009 . Genetic variability, shell and sperm morphology suggest that the surf clams Donax marincovichi and D. obesulus are one species . Journal of Molluscan Studies 75 : 381 – 390 . Google Scholar CrossRef Search ADS Coan EV , Scott PV , Bernard FR . 2000 . Bivalve seashells of western North America: marine bivalve mollusks from Arctic Alaska to Baja California. Santa Barbara Museum of Natural History Monographs Number 2; Studies in Biodiversity Number 2 . Santa Barbara : Santa Barbara Museum of Natural History . Congleton WR , Vassiliev T , Bayer RC , Pearce BR , Jacques J , Gillman C . 2006 . Trends in Maine softshell clam landings . Journal of Shellfish Research 25 : 475 – 480 . Google Scholar CrossRef Search ADS Crocetta F , Turolla E . 2011 . Mya arenaria Linné, 1758 (Mollusca: Bivalvia) in the Mediterranean Sea: its distribution revisited . Journal of Biological Research-Thessaloniki 16 : 188 – 193 . Cross ME , Bradley CR , Cross TF , Culloty S , Lynch S , McGinnity P , O’Riordan RM , Vartia S , Prodohl PA . 2016 . Genetic evidence supports recolonisation by Mya arenaria of western Europe from North America . Marine Ecology Progress Series 549 : 99 – 112 . Google Scholar CrossRef Search ADS Darriba D , Taboada GL , Doallo R , Posada D . 2012 . jModelTest 2: more models, new heuristics and parallel computing . Nature Methods 9 : 772 . Google Scholar CrossRef Search ADS PubMed Drozdov AL , Sharina SN , Tyurin SA . 2009 . Sperm ultrastructure in representatives of six bivalve families from Peter the Great Bay, Sea of Japan . Russian Journal of Marine Biology 35 : 236 – 241 . Google Scholar CrossRef Search ADS Dzyuba SM , Maslennikova LA . 1987 . Reproductive-cycle of the bivalve mollusk Mya japonica in Peter the Great Bay of the Sea of Japan . Biologiya Morya-Marine Biology 2 : 38 – 41 . Eckelbarger KJ , Young CM . 1999 . Ultrastructure of gametogenesis in a chemosynthetic mytilid bivalve (Bathymodiolus childressi) from a bathyal, methane seep environment (northern Gulf of Mexico) . Marine Biology 135 : 635 – 646 . Google Scholar CrossRef Search ADS Edgecombe GD , Giribet G , Wheeler WC . 2002 . Phylogeny of Henicopidae (Chilopoda: Lithobiomorpha): a combined analysis of morphology and five molecular loci . Systematic Entomology 27 : 31 – 64 . Google Scholar CrossRef Search ADS Elderkin CL , Clewing C , Ndeo OW , Albrecht C . 2016 . Molecular phylogeny and DNA barcoding confirm cryptic species in the African freshwater oyster Etheria elliptica Lamarck, 1807 (Bivalvia: Etheriidae) . Biological Journal of the Linnean Society 118 : 369 – 381 . Google Scholar CrossRef Search ADS Fofonoff PW , Ruiz GM , Steves B , Carlton JT . 2003 . California Non-native Estuarine and Marine Organisms (Cal-NEMO) system . Available at: http://invasions.si.edu/nemesis (accessed 2 May 2017 ). Folmer O , Black M , Hoeh W , Lutz R , Vrijenhoek R . 1994 . DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates . Molecular Marine Biology and Biotechnology 3 : 294 – 299 . Google Scholar PubMed Fujie T . 1957 . On the myarian pelecypods of Japan: part I. Summary of the study of the genus Mya from Hokkaido . Journal of the Faculty of Science, Hokkaido University. Series 4, Geology and Mineralogy 9 : 381 – 413 . Fujie T . 1962 . On the myarian pelecypods of Japan: Part II. Geological and geographical distribution of fossil and recent species, genus Mya . Journal of the Faculty of Science, Hokkaido University. Series 4, Geology and Mineralogy 11 : 399 – 429 . Geraghty J , Russell MP , Dollahon N . 2008 . A quantitative assessment of spermatozoan morphology in Nutricola confusa and Nutricola tantilla (Bivalvia: Veneridae) . Veliger 50 : 263 – 268 . Gomoiu M-T , Alexandrov B , Shadrin N , Zaitsev Y . 2002 . The Black Sea – a recipient, donor and transit area for alien species . In: Leppäkoski E , Gollasch S , Olenin S , eds. Invasive aquatic species of Europe. Distribution, impacts and management . Dordrecht : Kluwer Academic Publishers , 341 – 350 . Google Scholar CrossRef Search ADS Guindon S , Dufayard JF , Lefort V , Anisimova M , Hordijk W , Gascuel O . 2010 . New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0 . Systematic Biology 59 : 307 – 321 . Google Scholar CrossRef Search ADS PubMed Gwo JC , Liou CH , Cheng CH . 1996 . Ultrastructure of the spermatozoa of the Pacific oyster Crassostrea gigas (Mollusca, Bivalvia, Ostreidae) . Journal of Submicroscopic Cytology and Pathology 28 : 395 – 400 . Habe T . 1955 . Fauna of Akkeshi Bay. XXI. Pelecypoda and Scaphopoda . Publications of the Akkeshi Marine Biological Station 4 : 1 – 31 . Habe T . 1977 . Systematics of Mollusca in Japan. Bivalvia and Scaphopoda . Tokyo : Hokuryukan . Hillis DM , Bull JJ . 1993 . An empirical test of bootstrapping as a method for assessing confidence in phylogenetic analysis . Systematic Biology 42 : 182 – 192 . Google Scholar CrossRef Search ADS Hodgson AN , Bernard RTF , Vanderhorst G . 1990 . Comparative spermatology of three species of Donax (Bivalvia) from South-Africa . Journal of Molluscan Studies 56 : 257 – 265 . Google Scholar CrossRef Search ADS Hodgson AN , Devilliers CJ , Bernard RTF . 1987 . Comparative spermatology of two morphologically similar species of Solen (Mollusca, Bivalvia) . South African Journal of Zoology 22 : 264 – 268 . Google Scholar CrossRef Search ADS Huber M . 2010 . Compendium of bivalves. A full-color guide to 3’300 of the world’s marine bivalves: a status on Bivalvia after 250 years of research . Hackenheim : ConchBooks . Introini GO , Passos FD , Recco-Pimentel SM . 2013 . Comparative study of sperm ultrastructure of Donax hanleyanus and Donax gemmula (Bivalvia: Donacidae) . Acta Zoologica 94 : 261 – 266 . Google Scholar CrossRef Search ADS ITIS . 2017 . Mya Linnaeus, 1758 . Integrated Taxonomic Information System . Available at: https://www.itis.gov/servlet/SingleRpt/SingleRpt?search_topic=TSN&search_value=81691#null (accessed 10 May 2017 ). Jay JC . 1856 . Report on the shells collected by the Japan Expedition together with a list of Japan shells . In: Perry MC , Hawks FL , eds. Narrative of the expedition of an American squadron to the China Sea and Japan, performed in the year 1852, 1853 and 1854, under the command of Commodore M. C. Perry, United States Navy, by order of the Government of the United States . Washington : Tucker , 289 – 297 . Katoh K , Standley DM . 2013 . MAFFT multiple sequence alignment software version 7: improvements in performance and usability . Molecular Biology and Evolution 30 : 772 – 780 . Google Scholar CrossRef Search ADS PubMed Kienberger K , Carmona L , Pola M , Padula V , Gosliner TM , Cervera JL . 2016 . Aeolidia papillosa (Linnaeus, 1761) (Mollusca: Heterobranchia: Nudibranchia), single species or a cryptic species complex? a morphological and molecular study . Zoological Journal of the Linnean Society 177 : 481 – 506 . Google Scholar CrossRef Search ADS Kim JH , Chung JS , Park YJ . 2011 . Ultrastructures of germ cells during spermatogenesis and taxonomic values in sperm morphology in male Mya arenaria oonogai (Heterodonta: Myidae) . Korean Journal of Malacology 27 : 377 – 386 . Google Scholar CrossRef Search ADS Kwon OK , Min DK , Lee JR , Lee JS , Je JG , Choe BL . 2001 . Korean mollusks with color illustration . Seoul : Hanguel Publishing Company . Kwon OK , Park GM , Lee JS . 1993 . Colored shells of Korea . Seoul : Academy Publishing Co . Lasota R , Pierscieniak K , Garcia P , Simon-Bouhet B , Wolowicz M . 2016 . Large-scale mitochondrial COI gene sequence variability reflects the complex colonization history of the invasive soft-shell clam, Mya arenaria (L.) (Bivalvia) . Estuarine, Coastal and Shelf Science 181 : 256 – 265 . Google Scholar CrossRef Search ADS Laursen D . 1966 . The genus Mya in the Arctic Region . Malacologia 3 : 399 – 418 . Lavoie R . 2005 . Oyster culture in North America history, present and future . The 1st International Oyster Symposium Proceedings 19 : 9 . Layton KK , Martel AL , Hebert PD . 2014 . Patterns of DNA barcode variation in Canadian marine molluscs . PLoS ONE 9 : e95003 . Google Scholar CrossRef Search ADS PubMed Lutaenko KA , Noseworthy RG . 2012 . Catalogue of the living Bivalvia of the continental coast of the Sea of Japan (East Sea) . Vladivostok : Dalnauka . MacNeil FS . 1965 . Evolution and distribution of the genus Mya, and Tertiary migrations of Mollusca . United States, Department of the Interior, Geological Survey, Professional Paper 483G : 1 – 51 . Makiyama J . 1934 . The Asagaian molluscs of Yotukara and Matchgar . Memoirs of the College of Science Kyoto Imperial University, Series B 10 : 156 – 159 . Makiyama J . 1935 . The fossils of the genus Mya . Warera no Kobutsu 4 : 135 – 139 . Marko PB . 2002 . Fossil calibration of molecular clocks and the divergence times of geminate species pairs separated by the Isthmus of Panama . Molecular Biology and Evolution 19 : 2005 – 2021 . Google Scholar CrossRef Search ADS PubMed Nagao T , Inoue T . 1941 . Myarian fossils from the Cenozoic deposits of hokkaidô and Karahuto . Journal of the Faculty of Science, Hokkaido Imperial University. Series 4, Geology and Mineralogy 6 : 143 – 158 . Nakashima R . 1999 . Mya (Bivalvia: Myidae) from the Upper Miocene and Pliocene formations in Hokkaido, Northern Japan . Venus: Japanese Journal of Malacology 58 : 201 – 216 . Navarro D . 2014 . Learning statistics with R: a tutorial for psychology students and other beginners, Version 0.5 . Adelaide : University of Adelaide . Nikula R , Strelkov P , Väinölä R . 2007 . Diversity and trans-arctic invasion history of mitochondrial lineages in the North Atlantic Macoma balthica complex (Bivalvia: Tellinidae) . Evolution 61 : 928 – 941 . Google Scholar CrossRef Search ADS PubMed Okutani T . 2000 . Marine Mollusks in Japan . Tokyo : Tokai University Press . Padial JM , Miralles A , De la Riva I , Vences M . 2010 . The integrative future of taxonomy . Frontiers in Zoology 7 : 16 . Google Scholar CrossRef Search ADS PubMed Petersen GH . 1999 . Five recent Mya species, including three new species and their fossil connections . Polar Biology 22 : 322 – 328 . Google Scholar CrossRef Search ADS Pfenninger M , Reinhardt F , Streit B . 2002 . Evidence for cryptic hybridization between different evolutionary lineages of the invasive clam genus Corbicula (Veneroida, Bivalvia) . Journal of Evolutionary Biology 15 : 818 – 829 . Google Scholar CrossRef Search ADS Porter RG . 1974 . Reproductive cycle of the soft-shell clam, Mya arenaria, at Skagit Bay, Washington . Fishery Bulletin 72 : 648 – 656 . Powers SP , Bishop MA , Grabowski JH , Peterson CH . 2006 . Distribution of the invasive bivalve Mya arenaria L. on intertidal flats of southcentral Alaska . Journal of Sea Research 55 : 207 – 216 . Google Scholar CrossRef Search ADS Puillandre N , Lambert A , Brouillet S , Achaz G . 2012 . ABGD, Automatic Barcode Gap Discovery for primary species delimitation . Molecular Ecology 21 : 1864 – 1877 . Google Scholar CrossRef Search ADS PubMed R Core Team . 2015 . R: a language and environment for statistical computing . Vienna : R Foundation for Statistical Computing . Ratnasingham S , Hebert PD . 2007 . Bold: the barcode of life data system (http://www.barcodinglife.org) . Molecular Ecology Notes 7 : 355 – 364 . Google Scholar CrossRef Search ADS PubMed Rees DJ , Dioli M , Kirkendall LR . 2003 . Molecules and morphology: evidence for cryptic hybridization in African Hyalomma (Acari: Ixodidae) . Molecular Phylogenetics and Evolution 27 : 131 – 142 . Google Scholar CrossRef Search ADS PubMed Reynolds ES . 1963 . The use of lead citrate at high pH as an electron-opaque stain in electron microscopy . The Journal of Cell Biology 17 : 208 – 212 . Google Scholar CrossRef Search ADS PubMed Ronquist F , Teslenko M , van der Mark P , Ayres DL , Darling A , Höhna S , Larget B , Liu L , Suchard MA , Huelsenbeck JP . 2012 . MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space . Systematic Biology 61 : 539 – 542 . Google Scholar CrossRef Search ADS PubMed Scarlato OA . 1981 . Bivalve mollusks of temperate latitudes of the western portion of the Pacific Ocean. Guide books on the fauna of the USSR . Leningrad : Zoological Institute, USSR Academy of Sciences . Schlick-Steiner BC , Steiner FM , Seifert B , Stauffer C , Christian E , Crozier RH . 2010 . Integrative taxonomy: a multisource approach to exploring biodiversity . Annual Review of Entomology 55 : 421 – 438 . Google Scholar CrossRef Search ADS PubMed Schneider CA , Rasband WS , Eliceiri KW . 2012 . NIH Image to ImageJ: 25 years of image analysis . Nature Methods 9 : 671 – 675 . Google Scholar CrossRef Search ADS PubMed Strasser M . 1999 . Mya arenaria – an ancient invader of the North Sea coast . Helgoländer Meeresuntersuchungen 52 : 309 – 324 . Google Scholar CrossRef Search ADS Strasser CA , Barber PH . 2009 . Limited genetic variation and structure in softshell clams (Mya arenaria) across their native and introduced range . Conservation Genetics 10 : 803 – 814 . Google Scholar CrossRef Search ADS Strauch F . 1972 . Phylogenese, adaptation und migration einiger nordischer mariner Molluskengenera (Neptunea, Panomya, Cyrtodaria und Mya) . Abhandlungen der Senckenbergischen Naturforschenden Gesellschaft 531 : 1 – 211 . Streftaris N , Zenetos A . 2006 . Alien marine species in the Mediterranean – the 100 ‘worst invasives’ and their impact . Mediterranean Marine Science 7 : 87 – 118 . Google Scholar CrossRef Search ADS Tamura K , Stecher G , Peterson D , Filipski A , Kumar S . 2013 . MEGA6: molecular evolutionary genetics analysis version 6.0 . Molecular Biology and Evolution 30 : 2725 – 2729 . Google Scholar CrossRef Search ADS PubMed Tchang S , Qi Z , Li J . 1955 . Economic marine Mollusca in North China Sea . Beijing : Science Press . Torroglosa ME , Gimenez J . 2015 . Sperm ultrastructure in two species of Brachidontes (Bivalvia, Mytilidae) from the south-western Atlantic Ocean . Journal of the Marine Biological Association of the United Kingdom 95 : 991 – 998 . Google Scholar CrossRef Search ADS Vaidya G , Lohman DJ , Meier R . 2011 . SequenceMatrix: concatenation software for the fast assembly of multi-gene datasets with character set and codon information . Cladistics 27 : 171 – 180 . Google Scholar CrossRef Search ADS Väinölä R . 2003 . Repeated trans-Arctic invasions in littoral bivalves: molecular zoogeography of the Macoma balthica complex . Marine Biology 143 : 935 – 946 . Google Scholar CrossRef Search ADS Vermeij GJ . 1989 . Invasion and extinction: the last three million years of North Sea pelecypod history . Conservation Biology 3 : 274 – 281 . Google Scholar CrossRef Search ADS Whiting MF , Carpenter JC , Wheeler QD , Wheeler WC . 1997 . The Strepsiptera problem: phylogeny of the holometabolous insect orders inferred from 18S and 28S ribosomal DNA sequences and morphology . Systematic Biology 46 : 1 – 68 . Google Scholar PubMed Xiong B , Kocher TD . 1991 . Comparison of mitochondrial DNA sequences of seven morphospecies of black flies (Diptera: Simuliidae) . Genome 34 : 306 – 311 . Google Scholar CrossRef Search ADS PubMed Yoo JS . 1967 . Korean shells in colour . Seoul : Il Ji Sa Publishing Co . Zhang J , Xu F , Liu R . 2012 . The Myidae (Mollusca, Bivalvia) from Chinese waters with description of a new species . Zootaxa 3383 : 39 – 60 . Google Scholar CrossRef Search ADS © 2018 The Linnean Society of London, Zoological Journal of the Linnean Society

Journal

Zoological Journal of the Linnean SocietyOxford University Press

Published: Jan 17, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off