\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,greek,russian,english]{babel}



\usepackage{float}








\begin{document}

\title{Phylogeography of Mytilisepta virgata (Mytilidae: Bivalvia) in the
northwestern Pacific: cryptic mitochondrial lineages and mito-nuclear
discordance}



\author[1]{Yucheol Lee}%
\author[1]{Gang Ni}%
\author[1]{Jinkyung Shin}%
\author[1]{Taeho Kim}%
\author[1]{Elizabeth  Kern}%
\author[1]{Yuseob Kim}%
\author[2]{Seung-Chul Kim}%
\author[3]{Benny Chan}%
\author[4]{Ryutaro Goto}%
\author[5]{Tomoyuki Nakano}%
\author[1]{Joong-Ki Park}%
\affil[1]{Ewha Womans University}%
\affil[2]{Sungkyunkwan University College of Natural Science}%
\affil[3]{Academia Sinica}%
\affil[4]{Kyoto University}%
\affil[5]{Affiliation not available}%


\vspace{-1em}



  
  \date{May 21, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
The purplish bifurcate mussel Mytilisepta virgata is widely distributed
and represents one of the major components of the intertidal community
in the northwestern Pacific (NWP). Here, we characterized population
genetic structure of NWP populations throughout their whole distribution
range using both mitochondrial (mtDNA cox1) and nuclear (ITS1) markers.
Population genetic analyses for mtDNA cox 1 sequences revealed two
monophyletic lineages (i.e., southern and northern lineages)
geographically distributed according to the two different surface water
temperature zones in the NWP. The timing of the lineage split is
estimated at the Pliocene- mid-Pleistocene (5.49-1.61 Mya), which is
consistent with the timing of the historical isolation of the East
Sea/Sea of Japan from the South and East China Seas caused by sea level
decline during glacial cycles. Historical sea level fluctuation during
the Pliocene-Pleistocene and subsequent adaptation of mussels to
different surface water temperature zones may have contributed to
shaping the contemporary genetic diversity and deep divergence of the
two mitochondrial lineages. Unlike mtDNA sequences, a clear lineage
splitting between the two mitochondrial lineages was not found in ITS1
sequences, showing a star-like structure that is composed of a mixture
of southern and northern mitochondrial lineages. Possible scenarios are
proposed to explain this type of mito-nuclear discordance: stochastic
divergence in the coalescent processes of the two molecular markers, or
balancing selection under different marine environments. Future work is
required to address whether the thermal physiology of these mussels
correlates with the deep divergence of their mitochondrial genes.%
\end{abstract}%



\sloppy


\textbf{Introduction}

The northwestern Pacific (NWP) is a vast marine realm characterized by
unique tectonic and hydrologic features (Wang, 1999), providing a
natural setting in which to investigate mechanisms and processes leading
to marine diversification and speciation (Kong, Matsukuma, Hayashi,
Takada, \& Li, 2012; Ni, Li, Kong, \& Yu, 2014; Ho, Kwan, Kim, \& Won,
2015; Wang, Tsang, \& Dong, 2015; Qiu, Li, Lin, Ding, \& Miyamoto,
2016). Three continuous marginal seas are distributed from north to
south in the NWP, including the East Sea (ES)/Sea of Japan (SJ), East
China Sea (ECS) and South China Sea (SCS) (Tamaki \& Honza, 1991). The
size and connectivity of these seas probably dramatically changed
throughout the Pleistocene (Wang et al., 1997). During past glacial
epochs, sea levels dropped more than 120 m below the present level,
resulting in the exposure of sills and continental shelves in the NWP
(Wang \& Sun, 1994; Voris, 2000) and effectively closing the sea passage
(the Korean Strait, depth \textless{}130 m) between the ES/SJ and ECS,
resulting in their geographic separation (Kitamura, Takano, Takata, \&
Omote, 2001). The ECS experienced pronounced seawards migration
(\textasciitilde{}1200 km), and was reduced to an elongated and closed
trough (the Okinawa basin) (Xie, Jian, \& Zhao, 1995), while the SCS
became a semi-closed gulf with an area that was half its present size,
connected to the Pacific through the Bashi Strait (Wang, 1999).

The dramatic paleohydrogeological changes in the NWP probably had a
profound impact on marine species' distribution and evolution (Kong \&
Li, 2009; Liu, Li, Kong, \& Zheng, 2011). An evolutionary paradigm
emerging from marginal sea separation is that each sea basin served as
an independent refugium for surviving populations during the glacial
periods, and enhanced vicariant divergence among separated populations
(Liu, Gao, Wu, \& Zhang, 2007; Xu, Chan, Tsang, \& Chu, 2009; Ni et al.,
2014; Qiu et al., 2016). This pattern has been demonstrated in various
marine faunas, including fishes (Shen, Jamandre, Hsu, Tzeng, \& Durand,
2011; Song, Ma, Zhang, Gao, \& Sun, 2014; Qiu et al., 2016), molluscs
(Kong \& Li, 2009; Liu et al., 2011; Ni, Li, Kong, \& Zheng, 2012),
barnacles (Chan, Tsang \& Chu, 2007; Tsang et al. 2012) and mitten crab
(Xu et al., 2009; Xu \& Chu, 2012), all of which show substantial
intraspecific evolutionary partitions among the seas. However, the
paradigm is not applicable to all cases and is often complicated by
biotic and/or abiotic factors, such as surface water temperate
gradients, ocean currents, and diverse species-specific life-history
traits (Dong, Han, Ganmanee, \& Wang, 2015; He, Mukai, Chu, Ma, \&
Zhang, 2015; Ni, Li, Ni, Kong, \& Yu, 2015; Li et al., 2017; Ni et al.,
2020).

The surface water temperature gradient along the main coastline of the
NWP is substantial. This coastline spans a broad latitudinal range
roughly from 15 to 45\selectlanguage{ngerman}° N and represents a transition zone connecting the
tropical Indo-West Pacific and the cold North Pacific (Chen \& Jiao,
1997; Briggs \& Bowen, 2012). A biogeographic line, which begins at the
Yangtze Estuary in China and extends northeastwards to Jeju Island
(Korea) and southern Japan (Fig. 1), separates the North Pacific
Temperate Biotic Region (characterized by cold-temperate fauna) and
Indo-West Pacific Warm-water Biotic Region (characterized by
tropical/subtropical fauna) (Liu, 2013; Ni, Kern, Dong, Li, \& Park,
2017). Species with narrow thermal tolerances are restricted to one side
of the boundary because of the steep temperature gradient across the
Yangtze Estuary (Zhang, Qi, Zhang, \& Ma, 1963; Liu, 2013): the annual
sea surface temperature is \textgreater{}20°C on the south side of the
estuary, and rapidly decreases to \textless{}15°C on the north side
(Johnson \& Boyer, 2015). The surface water temperature gradient in this
region has been reported as an influential driver that accelerates
genetic differentiation (e.g., seagrass \emph{Zostera japonica} , Zhang,
Zhou, Xue, \& Liu, 2016) or cryptic speciation between neighboring
populations (e.g., fish \emph{Mugil cephalus} , Shen et al., 2011).
However, whether this temperature gradient is a general driver of
genetic differentiation in other marine species, especially coastal
invertebrate species, is still largely unknown.

Ocean currents can also impact species' genetic patterns in this region
(Dong et al., 2012; Han, Yanagimoto, Zhang, \& Gao, 2012; He et al.,
2015; Li et al., 2017). The main current systems in the NWP include the
cyclonic circulation of the Kuroshio Current (KC) and coastal currents
(Su \& Yuan, 2005; Wei, Yu, Ran, \& Zang, 2011). An extension of the
North Equatorial Current, the KC is a major current that dramatically
affects this region's environment and ecology (Kao, Wu, Hsin, \& Dai,
2006; Fujikura, Lindsay, Kitazato, Nishida, \& Shirayama, 2010). It
originates east of the Philippines and flows north-eastwards through
Taiwan and the Ryukyu Islands to the Pacific coast of southern mainland
Japan (Su, 2004). The KC and its branch currents transport large volumes
of water, salt, and nutrients from low latitude tropics to the northern
reaches of the marginal seas (Chen, 1997; Heath, Zenitani, Watanabe,
Kimura, \& Ishida, 1998). This region is also characterized by diverse
coastal currents (Su \& Yuan, 2005; Wei et al., 2011), such as the Subei
Coastal Current and the China Coastal Current along the coast of China,
the Korean Coastal Current around the Korean Peninsula, and the Oyashio
Current (a cold subarctic ocean current that originates in the Bering
Sea flowing southwest off the Kuril Islands and meet the KS east of
northern Japan). Oceanographic currents may potentially enhance
population connectivity, often producing different outcomes from sea
surface temperature gradient. Lack of genetic structuring among distant
populations can be maintained by continuous gene flow of current-driven
larval dispersal (Guo et al., 2015) and rafting events (Nikula, Fraser,
Spencer, \& Waters, 2010).

The purplish bifurcate mussel \emph{Mytilisepta virgata} (Wiegmann,
1837), a member of the family Mytilidae, is widely distributed in the
NWP including Korea, Japan, mainland China and Taiwan and forms large
mussel beds, constituting one of the major intertidal rocky shore
communities (Iwasaki, 1994; Liu \& Morton, 1994; Morton, 1995; Okutani,
2000; Min, Lee, Koh, \& Je, 2004). These mussel beds contribute to
marine biodiversity in the intertidal zone by providing habitat and
shelter for diverse marine invertebrate species (August, 1985; Seed,
1996). This species spawns twice a year, in spring (from February to
March) and autumn (from September to October) (Morton, 1995). The larval
duration of this species' planktonic stage is 4 days (Sasaki, 1984),
which can facilitate gene flow over a certain distance. These traits
(distribution range and larval ecology) make this species ideal for
testing the impacts and their relative strength of historical sea level
fluctuations, surface water temperature, and ocean currents on marine
phylogeography. In this study, based on comprehensive sampling
throughout the mussel's entire distribution in the NWP, we assessed
multiple lines of evidence (including mitochondrial and nuclear genes as
well as morphological data) to (1) test whether there exist
evolutionarily divergent linages in the NWP and (2) investigate driving
factors that may have prompted cryptic divergence.

\textbf{Materials and methods}

Sampling and DNA sequencing

For mitochondrial \emph{cox1} and nuclear ITS1 sequencing, a total of
680 \emph{M. virgata} individuals were collected from the intertidal
and/or subtidal zone at 34 localities on the seashores of China, Japan,
Korea, and Taiwan, which covers almost the entire distribution range
of\emph{M. virgata} (Fig. 1; Table 1). MtDNA \emph{cox1} and ITS1 from
four individuals of \emph{M. keenae} were also sequenced and used as
outgroup sequences in the phylogenetic analyses. All samples were
deposited as voucher specimens in the Marine Mollusk Resource Bank of
Korea (MMRBK: http://mmrbk.org/). Total genomic DNA was extracted from
the adductor muscle of each individual using a DNeasy tissue kit
following the manufacturer's instructions. A partial fragment of the
mtDNA \emph{cox1} gene was PCR-amplified using the universal primers
jgLCO1490 and jgHCO2198 (Geller, Meyer, Parker, \& Hawk, 2013) in total
reaction volumes of 50 \selectlanguage{greek}μ\selectlanguage{english}L, which included 2 \selectlanguage{greek}μ\selectlanguage{english}L of template DNA, 36.75 \selectlanguage{greek}μ\selectlanguage{english}L
of D.W., 5 \selectlanguage{greek}μ\selectlanguage{english}L of 10x Taq buffer, 1 \selectlanguage{greek}μ\selectlanguage{english}L of each primer (10 pmol), 4 \selectlanguage{greek}μ\selectlanguage{english}L of
0.2mM dNTP mixture, and 0.25 \selectlanguage{greek}μ\selectlanguage{english}L of \emph{Taq} polymerase (TaKaTa
Ex\emph{Taq} ). The PCR amplification conditions were as follows: 1
cycle of denaturation at 95 for 3 min followed by 35 cycles of
denaturation at 94 for 30 s, annealing at 45 for 30 s, elongation at
72 for 1 min, and a final extension at 72 for 10 min. To
cross-reference the mtDNA\emph{cox1} results, the nuclear ITS1 gene was
PCR-amplified and sequenced for 101 individuals from 23 populations,
which were subsampled so that mtDNA \emph{cox1} haplotype diversity was
well-represented in the ITS1 group (Table 1). The nuclear ITS1 gene
fragment was PCR-amplified using ITS1 (5\selectlanguage{ngerman}´-TCC GTA GGT GAA CCT GCG G-3´)
(White, Bruns, Lee, \& Taylor, 1990) and 5.8S-R (5´-TCG ATG AAG AAC GCA
GC-3´) (Vilgalys \& Hester, 1990). The PCR mixtures and conditions for
the ITS1 gene fragment were the same as for mtDNA \emph{cox1} , except
that the annealing temperature was 55\selectlanguage{english}. PCR products were isolated from
1\% agarose gels using a QIAquick gel extraction kit (QIAGEN Valencia,
USA) following standard protocols and sequenced directly using an ABI
PRISM 3700 DNA analyzer (Applied Biosystems, Foster City, USA).

Phylogenetic analyses and population structure

DNA sequences were aligned using Clustal W (Thompson, Higgins, \&
Gibson, 1994) with default options in Geneious Pro v8.1.9 (Biomatters,
Auckland, New Zealand). Haplotypes were identified with a Bayesian
coalescent-based program performed in DnaSP v5.10.1 (Librado \& Rozas,
2009). All mtDNA \emph{cox1} and ITS1 haplotypes were deposited in
GenBank (\emph{cox1} : MN728370-MN728538; ITS1: MN728333-MN728363).
Phylogenetic relationships among haplotypes were inferred using maximum
likelihood (ML) analysis and Bayesian inference (BI) implemented in
RAxML v8.2.9 (Stamatakis, 2014) and MrBayes v3.2.6 (Huelsenbeck \&
Ronquist, 2001) using nucleotide sequences of \emph{M. keenae} (GenBank
Accession nos. for \emph{cox1} : MN728366 -- MN728369; GenBank Accession
nos. for ITS1: MN728364 -- MN728365) as an outgroup. The best fit models
for the two genes were determined using the Akaike information criterion
(AIC) in jModelTest 2.1.10 (number of substitution schemes = 11;
Darriba, Taboada, Doallo, \& Posada, 2012): GTR+I for \emph{cox1} and
TVM+G (0.47) for ITS1. The ML analysis was performed using 1,000
bootstrap replicates. The BI analysis was conducted with two independent
runs, each with four simultaneous Monte Carlo Markov chain (MCMC) runs
for 1 x 10\textsuperscript{6} generations. Trees were sampled in every
100 generations and posterior probability was estimated after removing
an initial 25\% of generations (reaching a stationarity of the chains)
as burn-in. Haplotype networks were constructed using PopArt v1.7 (Leigh
\& Bryant, 2015) software via the TCS method (Clement, Posada, \&
Crandall, 2000).

Molecular diversity indices such as haplotype diversity (\emph{h} ),
nucleotide diversity (\selectlanguage{greek}\emph{π} \selectlanguage{english}), and mean number of pairwise
differences (\emph{k} ) of each population were estimated with ARLEQUIN
v3.5 (Excoffier \& Lischer, 2010). Pairwise genetic divergence between
populations was evaluated using the fixation index (\emph{F}
\textsubscript{st}) (Excoffier, Smouse, \& Quattro, 1992). A
hierarchical analysis of molecular variance (AMOVA) was conducted to
estimate population structure among geographical groups, which were
determined by a spatial analysis of molecular variance performed with
SAMOVA v2 (Dupanloup, Schneider, \& Excoffier, 2002). These analyses
were carried out with 1,000 permutations after sequential Bonferroni
adjustments using ARLEQUIN.

Morphological analysis

To confirm the molecular results (i.e., the presence of two distinct
lineages detected from the molecular analysis), we examined shell
morphometrics of \emph{M. virgata} populations. A total of 15 shell
characters (supplementary Fig. 1) were measured from 50 randomly
selected individuals from each lineage (100 individuals in total).
Measured shell characters were log-transformed, and principal components
analysis (PCA) was performed in R-studio to assess morphological
differences between the two lineages.

Divergence time and historical demography

The divergence time between the two \emph{M. virgata}
\emph{cox1}lineages was estimated using a strict clock model and the
HKY+I model in BEAST v2.4 (Drummond \& Bouckaert, 2015). The mtDNA
\emph{cox1}molecular clock for this species have not been empirically
assessed and we applied two different fossil-based divergence rates (a
minimum rate of 0.7\% and a maximum rate of 2.4\%/Mya) for mtDNA
\emph{cox1}estimated from two different molluscan taxa (0.7\%-1.2\%/Mya
{[}the ark clam; Marko, 2002{]}) and 2.4\%/Mya; {[}teguline gastropod;
Hellberg \& Vacquier, 1999{]}). To estimate the timing of divergence
between the two\emph{M. virgata} lineages, two independent runs with an
MCMC chain length of 2 x 10\textsuperscript{8} were performed with
sampling every 1,000 generations. Log and tree files were combined with
LogCombiner v2.4.8 (Rambaut \& Drummond, 2014a) with a burn-in of 25\%.
The effective sample size (ESS) values of all parameters above 1,000
were obtained in Tracer v1.6 (Rambaut, Suchard, Xie, \& Drummond, 2014).
Trees were used to build a maximum clade credibility tree in
TreeAnnotator v2.4.8 (Rambaut \& Drummond, 2014b) and visualized in
FigTree v1.4.3 (Rambaut, 2016). For testing historical demographic
events, Tajima's \emph{D} and Fu's \emph{F\textsubscript{S}} statistics
were first calculated to test for neutrality using ARLEQUIN. We then
used the pairwise mismatch distribution to test population expansion
using DnaSP.

\textbf{Results}

Phylogenetic analyses and population structure

The mtDNA \emph{cox1} gene fragment (661 bp in size) was sequenced for a
total of 680 individuals collected from 34 locations across the
northwestern Pacific (Fig. 1; Table 1). A total of 169 mtDNA
\emph{cox1}haplotypes were identified. Of these, 149 variable sites and
62 single-nucleotide polymorphic sites were detected, and 87 sites were
parsimoniously informative. The \emph{cox1} phylogenetic trees using ML
and BI methods showed identical tree topology: they showed two
monophyletic, divergent lineages (termed here the northern and southern
mitochondrial lineages, according to their geographic distribution) with
high support values (100\% BP and 1.00 BPP) (Fig. 2A). The haplotype
network displayed a dumbbell shape, with 48 steps (mutations) connecting
two haplogroups with a ``star genealogy'' (Fig. 2B). Pairwise sequence
divergence between the two lineages was unexpectedly very high, ranging
from 7.4\%-9.2\%, but there was noticeably low divergence within the
northern (0.5\%) and southern lineages (0.6\%), respectively.
Distribution of the two lineages clearly corresponds to their geographic
origins: the northern lineage includes all Korean populations and the
northern Japan populations (JRK, JAM, JOG, JHS, JIB, JSM, JTZ), with a
distribution range mainly in the cold surface water North Pacific
region. In contrast, the southern lineage is composed of all Chinese
populations, the Taiwan population and the southern Japan populations
(JBK, JAK, JTK, JHD, JKK), distributed northeastward along the tropical
Indo-West Pacific region. In addition, six populations (JSH {[}8{]}, JSS
{[}9{]}, JKS {[}11{]}, JOI {[}16{]}, JKG {[}17{]}, JNG {[}18{]}) on the
southeastern and western sides of the Japanese archipelago contained
both lineages (Fig. 1).

In addition to the \emph{cox1} data, 31 ITS1 genotypes were obtained
from 101 individuals that contained 29 polymorphic sites. Their pairwise
sequence divergence was very low (0.6\%) compared to mtDNA \emph{cox1} .
In contrast to the \emph{cox1} results, the ML and BI tree for the ITS1
data showed only a single clade where all individuals were mixed, with
no segregation between individuals representing the northern and
southern mitochondrial lineages (Fig. 3A). Network analysis of the ITS1
haplotype identified two major haplogroups that were only a single step
(mutation) away from each other, each forming a star-like genealogy that
comprised a mixture of northern and southern mitochondrial
\emph{cox1}members (Fig. 3B). Unlike the deep divergence detected
between the two mtDNA \emph{cox1} lineages, genotype distribution of
ITS1 did not correspond to geographic location.

The NWP populations of \emph{M. virgata} showed high haplotype diversity
with overall \emph{h} values of 0.874 for \emph{cox1} and 0.824 for ITS1
(Table 1). Nucleotide diversity (\selectlanguage{greek}\emph{π} \selectlanguage{english}) and nucleotide differences
(\emph{k} ) for \emph{cox1} were high: 0.0417 and 27.539, respectively.
In the six populations where the two divergent mtDNA lineages coexist,
nucleotide diversity (\selectlanguage{greek}\emph{π} \selectlanguage{english}) and nucleotide difference (\emph{k} )
ranged from 0.0192 to 0.0418 and from 12.699 to 27.613, respectively,
much higher than all other populations in this study (\selectlanguage{greek}\emph{π} \selectlanguage{english}=0.0008
to 0.0043, \emph{k} =0.544 to 2.828). However, ITS1 sequences showed
very low nucleotide diversity (\selectlanguage{greek}\emph{π=}\selectlanguage{english} 0.0026) and nucleotide
differences (\emph{k=} 0.909). In the SAMOVA test using the \emph{cox1}
gene,\emph{M. virgata} populations were subdivided into the two regional
groups consistent with the results from the phylogenetic analysis: one
group comprises populations of the cold-temperate zone that includes the
ES, Yellow Sea and northern Pacific coast of Japan (population nos. 1-9.
19-31), and the other group is composed of populations in the
subtropical-tropical zone represented by the ECS and SCS and the
southern Pacific coast of Japan (population nos. 10-18, 32-34). Genetic
differentiation (\emph{F} \textsubscript{ST}) between populations
belonging to different regions was substantial and statistically
significant, ranging from 0.416 to 0.984 (Supplementary Table 1).
However, \emph{F} \textsubscript{ST} values between populations within
the same regional group were less than 0.223, and most of the
comparisons were not statistically significant. The hierarchical AMOVA
analyses with locality groups inferred from the result of the SAMOVA
analysis indicated that 90.54\% of the total genetic variation was
contributed by `among-groups' variation
(\selectlanguage{russian}Ф\textsubscript{CT}=0.905,\emph{P} \textless{}0.001), whereas genetic
variation from `among populations within groups' and `within
populations' was 1.07\% (Ф\textsubscript{SC}=0.113, \emph{P}
\textless{}0.001) and 8.40\% (Ф\textsubscript{ST}=0.916, \emph{P}
\textless{}0.001), respectively. The AMOVA analysis excluding the 6
`mixed' populations (contained individuals representing both lineages)
shows that most of the genetic variation was observed from `among
groups' (97.33 \%, Ф\textsubscript{CT}=0.973, \emph{P} \textless{}0.001;
see Table 2 for more details).

Morphological analysis

Principle component analysis (PCA) of 15 morphological characters
(Supplementary Fig. 1) extracted from 100 randomly-chosen individuals (n
= 50 from the northern lineage, and n = 50 from the southern lineage)
were examined (Fig. 4). The first two components (PC1 and PC2) accounted
for 78.81\% of the total variation. Ordination of PC1 and PC2 scores of
100 individuals delimited by \emph{cox1} lineage indicated that the two
lineages mostly overlap along both axes. Therefore, we conclude that
these shell morphological characters do not differ between the two
mitochondrial lineages, despite their substantial genetic divergence.

Divergence time and historical demography

The divergence time between the northern and southern lineages estimated
using the two mtDNA \emph{cox1} molecular clock settings (0.7\% and
2.4\%/Mya) was 5.49--1.61 Mya, corresponding to the Pliocene to early-
to mid-Pleistocene, each coalescing at 0.94--0.30 Mya (the northern
lineage) and 0.99--0.29 Mya (the southern lineage), respectively (Fig.
5). Neutrality test showed significant negative (-) values of
Tajima's\emph{D} and Fu's \emph{F\textsubscript{S}} in both lineages
(Table 3). The sequential mismatch distributions for the two lineages
displayed a skewed unimodal peak close to the Y-axis (Fig. 6).

\textbf{Discussion}

Deep mitochondrial divergence

Recent molecular-based phylogeography studies in the NWP have uncovered
some general patterns: marine species often harbor genetically divergent
lineages in the SCS and ES (Cheang, Chu, \& Ang Jr, 2010
{[}macroalga:\emph{Sargassum hemiphyllum} {]}; Shen et al., 2011
{[}flathead grey mullet: \emph{Mugil cephalus} {]}; Ni et al., 2012
{[}venus clam:\emph{Cyclina sinensis} {]}; Cheng \& Sha, 2017
{[}Japanese mantis shrimp: \emph{Oratosquilla oratoria} {]}; Wang, Kong,
Chen, Matsukuma, \& Li, 2017 {[}spotted hard clam: \emph{Meretrix
petechialis} {]}). Here, using phylogenetic and network analyses of
mtDNA \emph{cox1} sequence data, we identified two distinct
mitochondrial lineages within \emph{M. virgata:} a northern lineage
found in the cold surface-water temperature zone in the northern NWP
(the Korean coast, the northern part of the ES/JS and the northern
Pacific coast of Japan), and a southern lineage occupying the warm
surface-water temperature zone in the NWP (ECS, SCS and the southern
Pacific coast of Japan). The fixation index (\emph{F}
\textsubscript{st}) and AMOVA analysis both support the conclusion that
these lineages are genetically distinct and display a different
biogeographic pattern except for some Japanese populations (JSH, JSS,
JKS, JOI, JKG, JNG; Fig. 1 and Table 1).

Sequence divergence between the two lineages is extremely high
(7.4-9.2\%), while within-lineage distance is low (a maximum of 0.6\%).
DNA barcoding analyses of other molluscan species using mtDNA\emph{cox1}
have reported a mean intraspecific variation less than 2\% (Meyer \&
Paulay, 2005; Zou, Li, \& Kong, 2012). However, deeper intraspecific
divergence (\textgreater{}2\%) has also been detected in some Atlantic
and Pacific molluscan species (gastropods and bivalve species); such
deep divergence has been attributed to population subdivision by a
biogeographic partition between geographic lineages (Layton, Martel, \&
Hebert, 2014). Very high levels of \emph{cox1}divergence have also been
reported in NWP bivalve species (Liu et al., 2011 {[}\textgreater{}4.3\%
in \emph{Atrina pectinata} {]}; Wang et al., 2017 {[}5.85\% in
\emph{Meretrix petechialis} {]}). Such divergence has been interpreted
as molecular evidence of a physical barrier reducing gene flow.
Historical isolation in marginal seas, caused by low sea level during
past glaciations, might be the primary driver for deep intraspecific
genetic splits in the NWP (Layton et al., 2014).

Discordance between mitochondrial and nuclear markers

The internal transcribed spacer (ITS1) of the ribosomal DNA has widely
been used as a molecular marker for phylogenetic analysis at species
level and for species identification of closely related molluscan taxa,
including mussel species (Santaclara et al., 2006; Wood, Apte, MacAvoy,
\& Gardner, 2007). Nevertheless, unlike the deep divergence we uncovered
in mitochondrial sequences, nuclear DNA sequences (ITS1) were not highly
divergent (a maximum pairwise sequence divergence of 0.6\%) and in the
phylogenetic tree, all ITS1 sequences clustered together, irrespective
of their geographic origin (Fig. 3A). The ITS1 network showed two
genotype groups (a single step away from each other), each forming a
star-like genealogy made of an admixture of northern and southern
mitochondrial \emph{cox1} members (Fig. 3B). This type of mito-nuclear
discordance is often found in biogeography studies and can be explained
by many possible hypotheses (reviewed in Toews and Brelsford 2012).
Further investigation will be needed to define the precise cause of
mito-nuclear discordance in this species, but we discussed two
hypotheses in the case of \emph{M. virgata} : incomplete lineage sorting
and selection on mitochondrial genes.

The first hypothesis invokes the stochasticity of coalescent processes
at the two independent loci, resulting in divergent patterns. Incomplete
lineage sorting is often mentioned in studies that uncover mito-nuclear
discordance because nuclear markers have a larger effective population
size (i.e. slower coalescence) than mitochondrial DNA (Moore, 1995;
Pesole, Gissi, De Chirico, \& Saccone, 1999). However, in the present
study the difference between the \emph{cox1} and ITS1 results are so
extreme (very deep divergence versus negligible divergence) that the
incomplete lineage sorting hypothesis is inadequate given the inferred
demographic history. Incomplete lineage sorting can only occur if two
subpopulations maintain large population sizes before they merge into a
common ancestral population; in contrast, our mtDNA \emph{cox1} and ITS1
genealogies both support very rapid and recent population growth
following a population bottleneck. This is indicated by the genealogies
with very short external branches only (a ``star genealogy;'' see Figs.
2A, 2B and Fig. 3), the significant negative values in Tajima's
\emph{D}and Fu's \emph{F\textsubscript{S}} test (Table 3), and the
unimodal peak skewed over the Y-axis in the pairwise mismatch
distribution analysis for the two mitochondrial lineages (Figs. 6A and
6B). Furthermore, the coalescence time between northern and southern
\emph{cox1} lineages is estimated to be more than five times older than
those of within-lineage coalescences (Fig. 5). Given a small ancestral
population(s) and/or a very old splits of an ancestral population,
incomplete lineage sorting at ITS1 or a long waiting time prior to the
coalescence of the mtDNA lineages are unlikely to be the reason for
mito-nuclear discordance.

Another hypothesis worth considering is that natural selection may have
produced long-term, balanced mtDNA polymorphism. Although it is a
long-held assumption that most of the sequence variation in mtDNA is
selectively neutral (Ballard \& Whitlock 2004; Meiklejohn, Montooth, \&
Rand, 2007), it is also well documented that balancing selection can
occur under spatial heterogeneity (Hedrick, 2006; Scott et al., 2011).
We note that the spatial distribution of the northern and
southern\emph{cox1} lineages strongly correlates with surface seawater
temperature (Fig. 1). Therefore, we suggest that these two
divergent\emph{cox1} lineages are a proxy for two mitochondrial
genotypes adapted to different surface water temperature zones, possibly
due to divergence in the energy-metabolic functions in mitochondria. To
test for divergent natural selection on the mt genome, we compared the
whole mitochondrial genome sequences of two individuals representing the
northern and southern mitochondrial lineages. Calculating the ratio of
nonsynonymous (\emph{dN} ) over synonymous (\emph{dS} ) substitutions
for the 12 mitochondrial protein-coding genes of the complete
mitochondrial genomes revealed that none of the 12 protein-coding genes
(including\emph{cox1} , the genetic marker used in this study) showed a
significant excess of dN/dS \textgreater{} 1 (Supplementary Table 2).
These results do not provide supporting evidence for the selection on
the mitochondrial genome; however, we cannot rule out the possibility
that the two divergent lineages might still be maintained by balancing
selection. The functional difference between two lineages might result
from a small number of amino acid substitutions (not readily detectable
in this case) or gene regulatory sequence differences. Additional
mitochondrial genome sequencing that can cover the genetic diversity for
the two mitochondrial lineages is further required to address whether
deep divergence in mitochondrial genes is related to adaptation to
different surface water temperature zones in the NWP.

Temperature-driven mitochondrial lineage splitting

Phylogeographic patterns of marine organisms in the NWP are shaped by
the interplay of historical and/or contemporary oceanography, and
present-day gene flow, which is often correlated with species-specific
life-history features such as larval traits (Cheang et al., 2012). The
timing of lineage splitting event between the two \emph{M.
virgata}lineages is estimated to be 5.49-1.61 Mya. Although there is
substantial deviation between estimates using the lowest and highest
calibration bounds of the molecular clock, the estimated divergence time
coincides with the glacial periods from the Pliocene to Pleistocene
epochs. Sea level fluctuation during this time affected population
fragmentation and demographic changes for many coastal species,
resulting in population structuring and/or lineage splitting in the NWP
(Shen et al., 2011; Ludt \& Rocha, 2014).Cyclic alternations of
glacial-interglacial periods dramatically affected the seafloor
topography and sea level dropped up to 120 m below its present level
during the glacial periods (i.e., early Pleiostocene), resulting in
historical isolation in the NW Pacific marginal seas: Pleiostocene
isolation of the ES/JS, ECS, and SCS (Wang, 1999; Voris, 2000). The
ES/JS was semi-enclosed and prevailed by cold surface water flows from
the northern sea during the Northern Hemisphere glaciation and
mid-Pleiostocene transition (Tada, 1994; Kitamura et al., 2001;
Matsuzaki, Itaki, Tada, \& Kamikuri, 2018). Meanwhile, the SCS was
separated from the ECS (Ota, 1998; see Fig. 7 for details) by a large
land bridge extending between eastern China and Taiwan Island (Kimura,
2000). Peng et al. (2019) found that there is a post-LGM expansion
of\emph{M. virgatus} in Zhejiang, China (belonging to the southern
lineage in the present study), based on population genetic analysis of
mtDNA\emph{cox1} and 16S rDNA sequences.

In addition to historical geography in the NWP, some environmental
factors (e.g., ocean current, surface water temperature and salinity)
might act in combination with stochastic processes such as local
selection and/or lineage sorting through time in shaping contemporary
population structure (Palumbi, 1994; Miglietta, Faucci, \& Santini,
2011). In our results, the geographic distribution of the two
mitochondrial lineages agrees with two biogeographic regions with
different surface water temperature zones: the northern lineage occupies
the North Pacific Temperate Biotic Region which is affected by the
cold-water Liman Current flowing southward from the Sea of Okhotsk into
in the East Sea (the Sea of Japan) and the Oyashio Current (also known
as the Kurile Current) that flows south and meets the Kuroshio Current
off the eastern shore of Japan. In contrast, the southern lineages are
occupied by Chinese populations, Taiwan population and the southern
Japan populations, all belonging to the Indo-West Pacific warm-water
Biotic Region which is affected by the warm, northeasterly flowing
Kuroshio Current. It is interesting to note that some Japanese
populations where the two lineages co-occur (JSH, JSS, JKS, JOI, JKG,
JNG; Table 1 and Fig. 1) agree well with the region where two currents
with different surface water temperatures meet: JOI (16), JKG (17), and
JNG (18) are found on the western side of Japanese archipelago where the
cold-water carrying Liman Current (LC) and the warm-water carrying
Tsushima Warm Current (TSWC) meet in the North Pacific Temperate Biotic
Region, whereas JSH (8), JSS (9), and JKS (11) are found in southeastern
side of Japan where the Kuroshio Current (KC) mixes with the
southeasterly flowing, cold-water carrying Oyashio Current (OC) in the
Indo-West Pacific warm-water Biotic Region.

Present-day oceanic currents are among the most influential factors that
contribute to continuous gene flow over a wide geographic range in many
marine invertebrates, including molluscan species (Dong et al., 2012;
Guo et al., 2015; Li et al., 2017). \emph{M. virgata} has a planktonic
larval stage for 4 days, allowing them to disperse considerable
distances via northeasterly flowing ocean currents (Sasaki, 1984). The
present-day geographic distribution of the northern and southern
lineages is assumed to be due to current-driven long-distance dispersal
of planktonic larva from the northern and southern glacial refugia. Deep
lineage splitting and geographic distribution of the two mitochondrial
lineages uncovered in this study suggest that along with ocean currents,
surface water temperature preference and/or differential selection
between the respective ancestral populations of the southern and
northern lineages have shaped contemporary phylogeographic patterns: The
ancestral population of the northern lineage would have been under
selection for adaptation to colder water, while the southern lineage
might have adapted to a warmer environment. This thermal adaptation
likely restricted gene flow between the two regional populations by
reducing migrant fitness, resulting in lineage splitting along a
temperature gradient (Teske, Von der Heyden, McQuaid, \& Barker, 2011;
Teske et al., 2019). Such a vicariant lineage splitting following
long-term adaptation to different temperature zones has been reported in
some other marine species including mantis shrimp (Cheng \& Sha, 2017),
mitten crab (Xu et al., 2009), barnacles (Tsang et al. 2012), and
flathead mullet (Shen et al., 2011). We cautiously suggest that cryptic
mitochondrial lineage splitting in \emph{M. virgata} might be the result
of an adaptive response (i.e., a balancing selection) to different
environments, most likely different surface water temperature zones. An
extensive genetic survey of nuclear genes, along with in-depth analysis
of balancing selection from mitochondrial genome sequences, would be
needed to clarify how these two mitochondrial lineages have adapted in
response to thermal differences in the NWP.

\textbf{Acknowledgements}

This work was supported by the National Research Foundation of Korea
(NRF) grant funded by the Korea government (MSIT) (Nos.
2017R1A2B4011597, 2020R1A2C2005393).

\textbf{References}

August, P. (1985). Islands of \emph{Mytilus} as a habitat for small
intertidal animals: effect of island size on community
structure.\emph{Marine Ecology Progress Series} , 25, 71-81. doi:
10.3354/meps025071

Ballard, J. W. O., \& Whitlock, M. C. (2004). The incomplete natural
history of mitochondria. \emph{Molecular Ecology} , 13(4), 729--744.
doi: 10.1046/j.1365-294X.2003.02063.x

Briggs, J. C., \& Bowen, B. W. (2012). A realignment of marine
biogeographic provinces with particular reference to fish
distributions.\emph{Journal of Biogeography} , 39(1), 12-30. doi:
10.1111/j.1365-2699.2011.02613.x

Chan, B. K. K., Tsang, L. M., \& Chu, K. H. (2007). Morphological and
genetic differentiation of the acorn barnacle \emph{Tetraclita
squamosa}(Crustacea, Cirripedia) in East Asia and description of a new
species of\emph{Tetraclita.} \emph{Zoologica Sripta} 36, 79-91.
doi:10.1111/j.1463-6409.2007.00260.x

Cheang, C. C., Chu, K. H., \& Ang Jr, P. O. (2010). Phylogeography of
the marine macroalga \emph{Sargassum hemiphyllum} (Phaeophyceae,
Heterokontophyta) in northwestern Pacific. \emph{Molecular Ecology} ,
19(14), 2933-2948. doi: 10.1111/j.1365-294X.2010.04685.x

Cheang, C. C., Tsang, L. M., Ng, W. C., Williams, G. A., Chu, K. H., \&
Chan, B. K. K. (2012). Phylogeography of the cold-water
barnacle\emph{Chthamalus challengeri} in the north-western Pacific:
effect of past population expansion and contemporary gene flow.
\emph{Journal of Biogeography} , 39(10), 1819-1835. doi:
10.1111/j.1365-2699.2012.02742.x

Chen, C. T. (1997). The Kuroshio intermediate water is the major source
of nutrients on the East China Sea continental shelf.\emph{Oceanographic
Literature Review} , 5(44), 531.

Chen, D. G., \& Jiao, Y. (1997). Comparative study of marine fishes and
their distributions in the China and Japan seas. \emph{Journal of Ocean
University of Qingdao} , 27(3), 305-312.

Cheng, J., \& Sha, Z. L. (2017). Cryptic diversity in the Japanese
mantis shrimp \emph{Oratosquilla oratoria} (Crustacea: Squillidae):
Allopatric diversification, secondary contact and
hybridization.\emph{Scientific Reports} , 7, 1927. doi:
10.1038/s41598-017-02059-7

Clement, M., Posada, D. C. K. A., \& Crandall, K. A. (2000). TCS: a
computer program to estimate gene genealogies. \emph{Molecular Ecology}
, 9(10), 1657-1659. doi: 10.1046/j.1365-294x.2000.01020.x

Darriba, D., Taboada, G. L., Doallo, R., \& Posada, D. (2012).
jModelTest 2: more models, new heuristics and parallel
computing.\emph{Nature Methods} , 9(8), 772. doi: 10.1038/nmeth.2109

Dong, Y. W., Han, G. D., Ganmanee, M., \& Wang, J. (2015). Latitudinal
variability of physiological responses to heat stress of the intertidal
limpet \emph{Cellana toreuma} along the Asian coast. \emph{Marine
Ecology Progress Series} , 529, 107-119. doi: 10.3354/meps11303

Dong, Y. W., Wang, H. S., Han, G. D., Ke, C. H., Zhan, X., Nakano, T.,
\& Williams, G. A. (2012). The impact of Yangtze river discharge, ocean
currents and historical events on the biogeographic pattern
of\emph{Cellana toreuma} along the China coast. \emph{PLoS ONE} , 7(4),
e36178. doi: 10.1371/journal.pone.0036178

Drummond, A. J., \& Bouckaert, R. R. (2015). \emph{Bayesian evolutionary
analysis with BEAST} . United Kingdom: Cambridge University Press.

Dupanloup, I., Schneider, S., \& Excoffier, L. (2002). A simulated
annealing approach to define the genetic structure of
populations.\emph{Molecular Ecology} , 11(12), 2571-2581. doi:
10.1046/j.1365-294X.2002.01650.x

Excoffier, L., \& Lischer, H. E. (2010). Arlequin suite ver 3.5: a new
series of programs to perform population genetics analyses under Linux
and Windows. \emph{Molecular Ecology Resources} , 10(3), 564-567. doi:
10.1111/j.1755-0998.2010.02847.x

Excoffier, L., Smouse, P. E., \& Quattro, J. M. (1992). Analysis of
molecular variance inferred from metric distances among DNA haplotypes:
application to human mitochondrial DNA restriction data.\emph{Genetics}
, 131(2), 479-491.

Fujikura, K., Lindsay, D., Kitazato, H., Nishida, S., \& Shirayama, Y.
(2010). Marine biodiversity in Japanese waters. \emph{PLoS ONE} , 5(8),
e11836. doi: 10.1371/journal.pone.0011836

Geller, J., Meyer, C., Parker, M., \& Hawk, H. (2013). Redesign of PCR
primers for mitochondrial cytochrome c oxidase subunit I for marine
invertebrates and application in all\selectlanguage{english}-taxa biotic surveys.\emph{Molecular
Ecology Resources} , 13(5), 851-861. doi: 10.1111/1755-0998.12138

Guo, X., Zhao, D., Jung, D., Li, Q., Kong, L. F., Ni, G., \ldots{} Park,
J. K. (2015). Phylogeography of the Rock Shell \emph{Thais
clavigera}(Mollusca): evidence for long-distance dispersal in the
northwestern Pacific. \emph{PloS ONE} , 10(7), e0129715. doi:
10.1371/journal.pone.0129715

Han, Z., Yanagimoto, T., Zhang, Y., \& Gao, T. (2012). Phylogeography
study of \emph{Ammodytes personatus} in Northwestern Pacific:
Pleistocene isolation, temperature and current conducted secondary
contact. \emph{PLoS ONE} , 7(5), e37425. doi:
10.1371/journal.pone.0037425

He, L., Mukai, T., Chu, K. H., Ma, Q., \& Zhang, J. (2015).
Biogeographical role of the Kuroshio Current in the amphibious
mudskipper \emph{Periophthalmus modestus} indicated by mitochondrial DNA
data. \emph{Scientific Reports} , 5, 15645. doi: 10.1038/srep15645

Heath, M., Zenitani, H., Watanabe, Y., Kimura, R., \& Ishida, M. (1998).
Modelling the dispersal of larval Japanese sardine, \emph{Sardinops
melanostictus} , by the Kuroshio Current in 1993 and
1994.\emph{Fisheries Oceanography} , 7(3-4), 335-346. doi:
10.1046/j.1365-2419.1998.00076.x

Hedrick, P. W. (2006). Genetic polymorphism in heterogeneous
environments: the age of genomics. \emph{Annual Review of Ecology,
Evolution, and Systematics} , 37, 67-93. doi:
10.1146/annurev.ecolsys.37.091305.110132

Hellberg, M. E., \& Vacquier, V. D. (1999). Rapid evolution of
fertilization selectivity and lysin cDNA sequences in teguline
gastropods. \emph{Molecular Biology and Evolution} , 16(6), 839-848.
doi: 10.1093/oxfordjournals.molbev.a026168

Ho, P. T., Kwan, Y. S., Kim, B., \& Won, Y. J. (2015). Postglacial range
shift and demographic expansion of the marine intertidal
snail\emph{Batillaria attramentaria} . \emph{Ecology and Evolution} ,
5(2), 419-435. doi: 10.1002/ece3.1374

Huelsenbeck, J. P., \& Ronquist, F. (2001). MRBAYES: Bayesian inference
of phylogenetic trees. \emph{Bioinformatics} , 17(8), 754-755. doi:
10.1093/bioinformatics/17.8.754

Iwasaki, K. (1994). Distribution and bed structure of the two intertidal
mussels, \emph{Septifer virgatus} (Wiegmann) and \emph{Hormomya
mutabilis} (Gould). \emph{Publications of the Seto Marine Biological
Laboratory} , 36(4), 223-247. doi: 10.5134/176237

Johnson, D. R., \& Boyer, T. P. (2015). \emph{East Asian Seas Regional
Climatology (Version 2)} . National Centers for Environmental
Information, NOAA, dataset. Available from
https://www.nodc.noaa.gov/OC5/regional\_climate/EASclimatology/

Kao, S. J., Wu, C. R., Hsin, Y. C., \& Dai, M. (2006). Effects of sea
level change on the upstream Kuroshio Current through the Okinawa
Trough. \emph{Geophysical Research Letters} . 33(16), L16604. doi:
10.1029/2006GL026822

Kimura, M. (2000). Paleogeography of the Ryukyu Islands. \emph{Tropics}
, 10(1), 5--24. doi: 10.3759/tropics.10.5

Kitamura, A., Takano, O., Takata, H., \& Omote, H. (2001). Late
Pliocene--early Pleistocene paleoceanographic evolution of thes Sea of
Japan. \emph{Palaeogeography, Palaeoclimatology, Palaeoecology} ,
172(1-2), 81-98. doi: 10.1016/S0031-0182(01)00272-3

Kong, L., Matsukuma, A., Hayashi, I., Takada, Y., \& Li, Q. (2012).
Taxonomy of \emph{Macridiscus} species (Bivalvia: Veneridae) from the
western Pacific: insight based on molecular evidence, with description
of a new species. \emph{Journal of Molluscan Studies} , 78(1), 1-11.
doi: 10.1093/mollus/eyr024

Kong, L. F., \& Li, Q. (2009). Genetic evidence for the existence of
cryptic species in an endangered clam \emph{Coelomactra antiquata}
.\emph{Marine Biology} , 156(7), 1507-1515. doi:
10.1007/s00227-009-1190-5

Layton, K. K., Martel, A. L., \& Hebert, P. D. (2014). Patterns of DNA
barcode variation in Canadian marine molluscs. \emph{PLoS ONE} , 9(4),
e95003. doi: 10.1371/journal.pone.0095003

Leigh, J. W., \& Bryant, D. (2015). Popart: full\selectlanguage{english}-feature software for
haplotype network construction. \emph{Methods in Ecology and Evolution}
, 6(9), 1110-1116. doi: 10.1111/2041-210x.12410

Li, J. J., Hu, Z. M., Gao, X., Sun, Z. M., Choi, H. G., Duan, D. L., \&
Endo, H. (2017). Oceanic currents drove population genetic connectivity
of the brown alga \emph{Sargassum thunbergii} in the north-west
Pacific.\emph{Journal of Biogeography} , 44(1), 230-242. doi:
10.1111/jbi.12856

Librado, P., \& Rozas, J. (2009). DnaSP v5: a software for comprehensive
analysis of DNA polymorphism data. \emph{Bioinformatics} , 25(11),
1451-1452. doi: 10.1093/bioinformatics/btp187

Liu, J. H., \& Morton, B. (1994). The temperature tolerances
of\emph{Tetraclita squamosa} (Crustacea: Cirripedia) and \emph{Septifer
virgatus} (Bivalvia: Mytilidae) on a sub\selectlanguage{english}-tropical rocky shore in Hong
Kong. \emph{Journal of Zoology} , 234(2), 325-339. doi:
10.1111/j.1469-7998.1994.tb06078.x

Liu, J., Li, Q., Kong, L. F., \& Zheng, X. D. (2011). Cryptic diversity
in the pen shell \emph{Atrina pectinata} (Bivalvia: Pinnidae): high
divergence and hybridization revealed by molecular and morphological
data. \emph{Molecular Ecology} , 20(20), 4332-4345. doi:
10.1111/j.1365-294X.2011.05275.x

Liu, J. X., Gao, T. X., Wu, S. F., \& Zhang, Y. P. (2007). Pleistocene
isolation in the Northwestern Pacific marginal seas and limited
dispersal in a marine fish, \emph{Chelon haematocheilus} (Temminck \&
Schlegel, 1845). \emph{Molecular Ecology} , 16(2), 275-288. doi:
10.1111/j.1365-294X.2006.03140.x

Liu, J. Y. (2013). Status of marine biodiversity of the China
Seas.\emph{PLoS ONE} , 8(1), e50719. doi: 10.1371/journal.pone.0050719

Ludt, W. B., \& Rocha, L. A. (2014). Shifting seas: the impacts of
Pleistocene sea-level fluctuations on the evolution of tropical marine
taxa. \emph{Journal of Biogeography} , 42(1), 25-38. doi:
10.1111/jbi.12416

Marko, P. B. (2002). Fossil calibration of molecular clocks and the
divergence times of geminate species pairs separated by the Isthmus of
Panama. \emph{Molecular Biology and Evolution} , 19(11), 2005-2021. doi:
10.1093/oxfordjournals.molbev.a004024

Matsuzaki, K. M., Itaki, T., Tada, R., \& Kamikuri, S. I. (2018).
Paleoceanographic history of the Japan Sea over the last 9.5 million
years inferred from radiolarian assemblages (IODP Expedition 346 Sites
U1425 and U1430). \emph{Progress in Earth and Planetary Science} , 5(1),
54. doi: 10.1186/s40645-018-0204-7

Meiklejohn, C. D., Montooth, K. L., \& Rand, D. M. (2007). Positive and
negative selection on the mitochondrial genome. \emph{Trends in
Genetics} . 23(6), 259--263. doi: 10.1016/j.tig.2007.03.008

Meyer, C. P., \& Paulay, G. (2005). DNA barcoding: error rates based on
comprehensive sampling. \emph{PLoS Biology} , 3(12), e422. doi:
10.1371/journal.pbio.0030422

Miglietta, M. P., Faucci, A., \& Santini, F. (2011). Speciation in the
sea: overview of the symposium and discussion of future
directions.\emph{Integrative and Comparative Biology} , 51(3), 449-455.
doi: 10.1093/icb/icr024

Min, D. K., Lee, J. S., Koh, D. B., \& Je, J. G. (2004). \emph{Mollusks
in Korea} . Seoul, Korea: Min Molluscan Research Institute.

Moore, W. S. (1995). Inferring phylogenies from mtDNA variation:
mitochondrial-gene trees versus nuclear-gene trees. \emph{Evolution} ,
49(4), 718-726. doi: 10.1111/j.1558-5646.1995.tb02308.x

Morton, B. (1995). The population dynamics and reproductive cycle
of\emph{Septifer virgatus} (Bivalvia: Mytilidae) on an exposed rocky
shore in Hong Kong. \emph{Journal of Zoology} , 235(3), 485-500. doi:
10.1111/j.1469-7998.1995.tb01764.x

Ni, G., Kern, E., Dong, Y. W., Li, Q., \& Park, J. K. (2017). More than
meets the eye: the barrier effect of the Yangtze River
outflow.\emph{Molecular Ecology} , 26(18), 4591-4602. doi:
10.1111/mec.14235

Ni, G., Kim, T., Shin, Y., Park, J., Lee, Y., Kil, H. J., \& Park, J. K.
(2020). Life-history features and oceanography drive phylogeographic
patterns of the chiton \emph{Acanthochitona} cf.
\emph{rubrolineata}(Lischke, 1873) in the northwestern Pacific.
\emph{PeerJ} , 8, e8794. doi: 10.7717/peerj.8794

Ni, G., Li, Q., Kong, L. F., \& Yu, H. (2014). Comparative
phylogeography in marginal seas of the northwestern
Pacific.\emph{Molecular Ecology} , 23(3), 534-548. doi:
10.1111/mec.12620

Ni, G., Li, Q., Kong, L. F., \& Zheng, X. D. (2012). Phylogeography of
bivalve \emph{Cyclina sinensis} : testing the historical glaciations and
Changjiang River outflow hypotheses in northwestern Pacific. \emph{PLoS
ONE} , 7(11), e49487. doi: 10.1371/journal.pone.0049487

Ni, G., Li, Q., Ni, L. H., Kong, L. F., \& Yu, H. (2015). Population
subdivision of the surf clam \emph{Mactra chinensis} in the East China
Sea: Changjiang River outflow is not the sole driver. \emph{PeerJ} , 3,
e1240. doi: 10.7717/peerj.1240

Nikula, R., Fraser, C., Spencer, H., \& Waters, J. (2010). Circumpolar
dispersal by rafting in two subantarctic kelp-dwelling
crustaceans.\emph{Marine Ecology Progress Series} , 405, 221-230. doi:
10.3354/meps08523

Okutani, T. (2000). \emph{Marine mollusks in Japan} . Tokyo, Japan:
Tokai University Press.

Ota, H. (1998). Geographic patterns of endemism and speciation in
amphibians and reptiles of the Ryukyu Archipelago, Japan, with special
reference to their paleogeographical implications. \emph{Researches on
Population Ecology} , 40(2), 189-204. doi: 10.1007/bf02763404

Palumbi, S. R. (1994). Genetic divergence, reproductive isolation, and
marine speciation. \emph{Annual Review of Ecology and Systematics} ,
25(1), 547-572. doi: 10.1146/annurev.es.25.110194.002555

Pesole, G., Gissi, C., De Chirico, A., \& Saccone, C. (1999). Nucleotide
substitution rate of mammalian mitochondrial genomes. \emph{Journal of
Molecular Evolution} , 48(4), 427-434. doi: 10.1007/pl00006487

Peng, X., Wu, X., Zhang, Y., Zhang, H., Tang, H., \& Kang, B. (2019)
Genetic diversity and metapopulation structures of two intertidal
species along the coast of Zhejiang, China: implications for
conservation, \emph{Mitochondrial DNA Part A} , 30,4, 674-681.
doi:10.1080/24701394.2019.1616702

Qiu, F., Li, H., Lin, H., Ding, S., \& Miyamoto, M. M. (2016).
Phylogeography of the inshore fish, \emph{Bostrychus sinensis} , along
the Pacific coastline of China. \emph{Molecular Phylogenetics and
Evolution} , 96, 112-117. doi: 10.1016/j.ympev.2015.11.020

Rambaut, A. (2016). \emph{FigTree v1.4.3} . Available from:
http://tree.bio.ed.ac.uk/software/figtree.

Rambaut, A., \& Drummond, A. J. (2014a). \emph{LogCombiner v.2.4.0} .
Part of the BEAST package. Available from: http://beast2.org.

Rambaut, A., \& Drummond, A. J., (2014b). \emph{TreeAnnotator v. 2.4.0}
. Part of the BEAST package. Available from: http://beast2.org.

Rambaut, A., Suchard, M. A., Xie, D., \& Drummond, A. J.
(2014).\emph{Tracer v1.6} . Available from:
http://beast.bio.ed.ac.uk/Tracer.

Santaclara, F. J., Espi\selectlanguage{ngerman}ñeira, M., Cabado, A. G., Aldasoro, A.,
Gonzalez-Lavín, N., \& Vieites, J. M. (2006). Development of a method
for the genetic identification of mussel species belonging
to\emph{Mytilus} , \emph{Perna} , \emph{Aulacomya} , and other
genera.\emph{Journal of Agricultural and Food Chemistry} , 54(22),
8461-8470. doi: 10.1021/jf061400u

Sasaki, R. (1984). Reproduction and larval development of \emph{Septifer
virgatas} {[}mussels{]} in Kesennuma bay, Miyagi prefecture
{[}Japan{]}.\emph{The Aquiculture} , 31, 214-219.

Scott G. R., Schulte, P. M., Egginton, S., Scott, A. L. M., Richards, J.
G., \& Milsom, W. K. (2011). Molecular evolution of cytochrome c oxidase
underlies high\selectlanguage{english}-altitude adaptation in the bar\selectlanguage{english}-headed
goose.\emph{Molecular Biology and Evolution} , 28(1), 351-363. doi:
10.1093/molbev/msq205

Seed, R. (1996). Patterns of biodiversity in the macro-invertebrate
fauna associated with mussel patches on rocky shores. \emph{Journal of
the Marine Biological Association of the United Kingdom} , 76(1),
203-210. doi: 10.1017/s0025315400029131

Shen, K. N., Jamandre, B. W., Hsu, C. C., Tzeng, W. N., \& Durand, J. D.
(2011). Plio-Pleistocene sea level and temperature fluctuations in the
northwestern Pacific promoted speciation in the globally-distributed
flathead mullet \emph{Mugil cephalus} . \emph{BMC Evolutionary Biology}
, 11(1), 83. doi: 10.1186/1471-2148-11-83

Song, N., Ma, G., Zhang, X., Gao, T., \& Sun, D. (2014). Genetic
structure and historical demography of \emph{Collichthys
lucidus}inferred from mtDNA sequence analysis. \emph{Environmental
Biology of Fishes} , 97(1), 69-77. doi: 10.1007/s10641-013-0124-8

Stamatakis, A. (2014). RAxML version 8: a tool for phylogenetic analysis
and post-analysis of large phylogenies. \emph{Bioinformatics} , 30(9),
1312-1313. doi: 10.1093/bioinformatics/btu033

Su, J. L. (2004). Overview of the South China Sea circulation and its
influence on the coastal physical oceanography outside the Pearl River
Estuary. \emph{Continental Shelf Research} , 24(16), 1745-1760. doi:
10.1016/j.csr.2004.06.005

Su, J. L., \& Yuan, Y. L. (2005). \emph{Coastal hydrology of China} .
Beijing, China: Ocean Press.

Tada, R. (1994). Paleoceanographic evolution of the Japan
Sea.\emph{Palaeogeography, Palaeoclimatology, Palaeoecology} , 108(3-4),
487-508. doi: 10.1016/0031-0182(94)90248-8

Tamaki, K., \& Honza, E. (1991). Global tectonics and formation of
marginal basins: role of the western Pacific. \emph{Episodes} , 14(3),
224-230. doi: 10.18814/epiiugs/1991/v14i3/005

Teske P. R., Sandoval-Castillo, J., Golla, T. R., Emami-Khoyi, A., Tine,
M., von der Heyden, S., \& Beheregaray, L. B. (2019). Thermal selection
as a driver of marine ecological speciation. \emph{Proceedings of the
Royal Society B} , 286(1869), 20182023. doi: 10.1098/rspb.2018.2023

Teske P. R., Von der Heyden, S., McQuaid C. D., \& Barker, N.P. (2011).
A review of marine phylogeography in southern Africa. \emph{South
African Journal of Science} , 107(5-6), 43-53. doi:
10.4102/sajs.v107i5/6.514

Thompson, J. D., Higgins, D. G., \& Gibson, T. J. (1994). CLUSTAL W:
improving the sensitivity of progressive multiple sequence alignment
through sequence weighting, position-specific gap penalties and weight
matrix choice. \emph{Nucleic Acids Research} , 22(22), 4673-4680. doi:
10.1093/nar/22.22.4673

Toews, D. P. L., \& Brelsford, A. (2012). The biogeography of
mitochondrial and nuclear discordance in animals. \emph{Molecular
Ecology} , 21(16), 3907--3930. doi: 10.1111/j.1365-294x.2012.05664.x

Tsang, L. M,, Wu, T. H., Shih, H. T., Williams, G. A., Chu, K. H., and
Chan, B. K. K. (2012) Genetic and morphological differentiation of the
Indo-West Pacific intertidal barnacle Chthamalus
malayensis.\emph{Integrative and Comparative Biology} , 52(3), 388--409.
doi:10.1093/icb/ics044

Vilgalys, R., \& Hester, M. (1990). Rapid genetic identification and
mapping of enzymatically amplified ribosomal DNA from
several\emph{Cryptococcus} species. \emph{Journal of Bacteriology} ,
172(8), 4238-4246. doi: 10.1128/jb.172.8.4238-4246.1990

Voris, H. K. (2000). Maps of Pleistocene sea levels in Southeast Asia:
shorelines, river systems and time durations. \emph{Journal of
Biogeography} , 27(5), 1153-1167. doi: 10.1046/j.1365-2699.2000.00489.x

Wang, J., Tsang, L. M., \& Dong, Y. W. (2015). Causations of
phylogeographic barrier of some rocky shore species along the Chinese
coastline. \emph{BMC Evolutionary Biology} , 15(1), 114. doi:
10.1186/s12862-015-0387-0

Wang, P. X. (1999). Response of western Pacific marginal seas to glacial
cycles: paleoceanographic and sedimentological features. \emph{Marine
Geology} , 156(1-4), 5-39. doi: 10.1016/s0025-3227(98)00172-8

Wang, P. X., Bradshaw, M., Ganzei, S. S., Tsukawaki, S., Binhassan, K.,
Hantoro, W. S., \ldots{} Kagami H. (1997). West Pacific marginal seas
during last glacial maximum: amplification of environmental signals and
its impact on monsoon climate. In \emph{Proceedings of the 30th
International Geological Congress} . Volume 13, pp. 65-86.

Wang, P. X., \& Sun, X. J. (1994). Last glacial maximum in China:
comparison between land and sea. \emph{Catena} , 23(3-4), 341-353. doi:
10.1016/0341-8162(94)90077-9

Wang, X., Kong, L., Chen, J., Matsukuma, A., \& Li, Q. (2017).
Phylogeography of bivalve \emph{Meretrix petechialis} in the
Northwestern Pacific indicated by mitochondrial and nuclear DNA
data.\emph{PloS ONE} , 12(8), e0183221. doi:
10.1371/journal.pone.0183221

Wei, Q., Yu, Z., Ran, X., \& Zang, J. (2011). Characteristics of the
western coastal current of the Yellow Sea and its impacts on material
transportation. \emph{Advances in Earth Science} , 26(2), 145-156.

White, T. J., Bruns, T., Lee, S., \& Taylor, J. L. (1990). Amplification
and direct sequencing of fungal ribosomal RNA genes for
phylogenetics.pdf. \emph{PCR protocols: a guide to methods and
applications,} 315-322. doi: 10.1016\%2Fb978-0-12-372180-8.50042-1

Wood, A. R., Apte, S., MacAvoy, E. S., \& Gardner, J. P. (2007). A
molecular phylogeny of the marine mussel genus \emph{Perna} (Bivalvia:
Mytilidae) based on nuclear (ITS1 and 2) and mitochondrial (COI) DNA
sequences. \emph{Molecular Phylogenetics and Evolution} , 44(2),
685-698. doi: 10.1016/j.ympev.2006.12.019

Xie, C., Jian, Z., \& Zhao, Q. (1995). Paleogeographic maps of the China
Seas at the last glacial maximum. In \emph{WESTPAC Paleogeographic
Maps}(p. 75). Shanghai, China: UNESCO/IOC Publications.

Xu, J., Chan, T. Y., Tsang, L. M., \& Chu, K. H. (2009). Phylogeography
of the mitten crab \emph{Eriocheir sensu stricto} in East Asia:
Pleistocene isolation, population expansion and secondary
contact.\emph{Molecular Phylogenetics and Evolution} , 52(1), 45-56.
doi: 10.1016/j.ympev.2009.02.007

Xu, J., \& Chu, K. H. (2012). Genome scan of the mitten
crab\emph{Eriocheir sensu stricto} in East Asia: Population
differentiation, hybridization and adaptive speciation. \emph{Molecular
Phylogenetics and Evolution} , 64(1), 118-129. doi:
10.1016/j.ympev.2012.03.009

Zhang, X., Qi, Z. Y., Zhang, F. S., \& Ma, X. T. (1963). A preliminary
study of the demarcation of marine molluscan faunal regions of China and
its adjacent waters. \emph{Oceanologia et Limnologia Sinica} , 2,
124-138.

Zhang, X., Zhou, Y., Xue, D. X., \& Liu, J. X. (2016). Genetic
divergence of the endangered seagrass \emph{Zostera japonica} Ascherson
\& Graebner between temperate and subtropical coasts of China based on
partial sequences of matK and ITS. \emph{Biochemical Systematics and
Ecology} , 68, 51-57. doi: 10.1016/j.bse.2016.06.008

Zou, S., Li, Q., \& Kong, L. (2012). Multigene barcoding and phylogeny
of geographically widespread muricids (Gastropoda: Neogastropoda) along
the coast of China. \emph{Marine Biotechnology} , 14(1), 21-34. doi:
10.1007/s10126-011-9384-5

\textbf{Data accesibility}

Mitochondrial \emph{cox1} and nuclear ITS1 sequences are available on
GenBank (accession nos. MN728370-MN728538 and MN728333-MN728363).

\textbf{Author contributions}

Study conception and design: JKP. Data collection: YL, JS, TK, SC, BC,
RG, YJW, YK, SCK, and TN. Data analysis and interpretation: YL, GN, YK,
and JKP. Manuscript preparation: YL, GN, EK, YK, and JKP.

Figure legends

\textbf{Figure 1.} Sample locations of \emph{M. virgata} and
distribution of the two mtDNA \emph{cox1} lineages. The map shows sea
surface temperature in the northwestern Pacific in spring (2003 --
2011). CCC, China Coastal Current; TWC, Taiwan Warm Current; CRDW,
Changjian River Diluted Water; SCC, Subei Coastal Current; YSWC, Yellow
Sea Warm Current; EKWC, East Korea Warm Current; TSWC, Tsushima Warm
Current; OC, Oyashio Current; LC, Liman Current. The dash/dotted line
separates two biogeographic regions: the North Pacific Temperate Biotic
Region; above the line) and Indo-West Pacific Warm-water Biotic Region;
below the line). Populations are labelled with numbers (the southern and
northern lineages are indicated in red and blue colored circles,
respectively) that correspond with those shown in Table 1.

\textbf{Figure 2.} Phylogenetic tree and haplotype network constructed
from mtDNA \emph{cox1} . In the phylogenetic tree (A), numbers on
branches are statistical support values for ML (bootstrap values)/BI
(posterior probabilities). In the haplotype network (B), colors
represent sampling locations of \emph{M. virgata} . The size of each
circle is proportional to the frequency of that haplotype. The crossbars
on the lines indicate the number of nucleotide substitutions. Open
circles represent hypothetical haplotypes.

\textbf{Figure 3.} Phylogenetic tree and haplotype network constructed
from nrDNA ITS1. In the phylogenetic tree (A), numbers at nodes are
statistical support values for ML (bootstrap values)/BI (posterior
probabilities). In the haplotype network (B), each circle represents one
allele and color indicates lineage population. Circle size is
proportional to the frequency of that haplotype. The crossbars on each
branch indicate the number of nucleotide substitution.

\textbf{Figure 4.} Plot of principal component analysis of northern and
southern lineages using 15 shell morphological characters. PC1 explains
72.14\% of the variation; PC2 6.67\% of the variation.

\textbf{Figure 5.} Estimates of divergence times using BEAST. A,
0.7\%/Mya; B, 2.4\%/Mya. Numbers above nodes are the mean divergence
times. Horizontal bars represent 95\% highest posterior densities of the
divergence time estimates.

\textbf{Figure 6.} Mismatch distributions for the northern and southern
lineages using mtDNA \emph{cox1} sequence data. The red dotted lines
represent the observed frequency of pairwise differences, and the blue
solid lines indicate the expected values of the sudden population
expansion model.

\textbf{Figure 7.} Hypothesized NW Pacific coastline (dotted line)
during early Pleistocene showing that three NWP marginal seas (the East
Sea {[}ES{]}/the Sea of Japan {[}SJ{]}, East China Sea {[}ECS{]} and
South China Sea {[}SCS{]}) were isolated due to sea level decline
(modified from Ota, 1998). The current land mass topography is indicated
in dark grey.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-1/Fig-1}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-2/Fig-2}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-3/Fig-3}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-4/Fig-4}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-5/Fig-5}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-6/Fig-6}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-7/Fig-7}

\end{center}
\end{figure}

\textbf{Hosted file}

\verb`Tables 1_2_3.docx`  available at \url{https://authorea.com/users/351687/articles/476224-phylogeography-of-mytilisepta-virgata-mytilidae-bivalvia-in-the-northwestern-pacific-cryptic-mitochondrial-lineages-and-mito-nuclear-discordance}

\selectlanguage{english}
\FloatBarrier
\end{document}

