100 
Fishery Bulletin 119(2—3) 
Cycler Dice Gradient (TaKaRa Bio, Shiga, Japan), each 
with 10 pL of reaction solution containing 10 ng template 
DNA, 0.5 pM of each primer, 0.2 mM dNTP mixture, 
1x PCR buffer (10 mM Tris—HCl pH 8.3, 50 mM KCl, 
and 1.5 mM MgCl,), and 0.5 U Ex Taq DNA polymerase 
(TaKaRa Bio). 
The reaction conditions were as follows: an initial 
denaturation at 94°C for 2 min; 30 cycles of denaturation 
at 94°C for 30 s, annealing at 62°C for 30 s, and exten- 
sion at 72°C for 1 min; and a final extension of 72°C for 
10 min. The PCR products, purified by using illustra Exo- 
Star (Cytiva, Malborough, MA), were sequenced on an 
Applied Biosystems ABI Prism 3130 XL Genetic Analyzer 
(Thermo Fisher Scientific Inc., Waltham, MA), with the 
Applied Biosystems BigDye Terminator Cycle Sequenc- 
ing Kit (vers. 3.1; Thermo Fisher Scientific Inc.) and the 
primers used in the PCR. Forward and reverse sequences 
were connected by using GENETYX, vers. 8 (GHNETYX 
Corp., Tokyo, Japan). All sequences have been deposited 
in the DNA Data Bank of Japan under accession numbers 
LC547089-LC547120. 
The DNA sequences and the published sequences of 
polkadot skate from waters around Jindo Island, off the 
southwestern coast of the Korean Peninsula, near the 
boundary of the East China Sea and the Yellow Sea (regis- 
tered as D. kwangtungensis; accession number KF318309; 
Jeong et al., 2015) and from Taiwan (the specific locality 
is unknown; accession number EU870495; Su and Hwang, 
available from the DNA Data Bank of Japan at website) 
were aligned with CLUSTAL W, vers. 2.1 (Larkin et al., 
2007). Haplotype diversity (h) (Nei, 1987) and nucleotide 
diversity (m) (Nei, 1987) were calculated for each geo- 
graphic population, meaning the population at each sam- 
pling location, by using Arlequin, vers. 3.5 (Excoffier and 
Lischer, 2010). 
Using the bigtail skate (LC547121; this study), the 
acutenose skate (accession number LC547122; this 
study), and the roughskin skate (D. trachydermus) (acces- 
sion number KR1526438; Vargas-Caro et al., 2016) as out- 
groups, we performed phylogenetic analysis with MEGA, 
vers. 6.0 (Tamura et al., 2013). The best nucleotide sub- 
stitution model was determined by using Bayesian infor- 
mation criterion. Subsequently, the maximum likelihood 
tree of haplotypes was constructed from evolutionary 
distances calculated by using the Tamura—Nei model 
(Tamura and Nei, 1993) with the G parameter (ensuring 
heterogeneous rates across sites following the gamma 
distribution) (TN93+G model). Confidence values were 
estimated through bootstrapping analysis (Felsenstein, 
1985) with 1000 replications. The time since lineage 
divergence (t) was estimated from net average distances 
between lineages (d,) (Nei and Li, 1979; Nei and Miller, 
1990) on the basis of the following relationship: t=d,/2p, 
where jp is the lineage mutation rate. Values of d, were 
calculated by using p-distances, the proportion (p) of 
nucleotide sites at which 2 sequences being compared are 
different, determined with MEGA. We assumed a diver- 
gence rate of 1.0—1.6% per million years, as this rate has 
been suggested for the mt cyt 6 gene of a species of Raja 
(Chevolot et al., 2006b). A minimum spanning network of 
haplotypes was constructed on the basis of the minimum 
sequence difference by using TCS, vers. 1.21 (Clement 
et al., 2000). 
A hierarchical analysis of molecular variance (AMOVA) 
(Excoffier et al., 1992) was run in Arlequin to estimate 
the partitioning of genetic variation within and between 
populations (or groups). We assumed different groupings 
of locations to find the optimal grouping of locations with 
the highest value of molecular genetic diversity among 
geographic groups (®c7). The pairwise values of molecu- 
lar genetic diversity among populations (®,7) were cal- 
culated, and their significance was tested with 10,000 
permutations and adjusted with a sequential Bonferroni 
correction (Rice, 1989). 
Analysis of mismatch distribution was performed with 
Arlequin to infer the recent demographic history of the 
samples and to test the fit of the observed distribution 
to that expected under the sudden expansion model. The 
goodness of fit between the observed and expected dis- 
tributions was tested by using the sum of squared devi- 
ations (Schneider and Excoffier, 1999) and Harpending’s 
raggedness index (Harpending, 1994). The demographic 
parameters of time since expansion in units of muta- 
tional time and population size before (85) and after (8,) 
a rapid change in population size in units of mutational 
time (Rogers and Harpending, 1992; Schneider and Excof- 
fier, 1999) were obtained by using a parametric boot- 
strap approach with 10,000 replications. Neutrality of 
the sequence variation was verified by using Tajima’s D 
statistic (Tajima, 1989a, 1989b) and Fu’s Fz test (Fu, 1997), 
conducted in Arlequin with 10,000 simulated samples. 
We constructed a Bayesian skyline plot (BSP) (Drum- 
mond et al., 2005) with the TN93+G model to estimate 
past dynamics in the female effective population size (NV,) 
multiplied by the generation interval through time, using 
BEAST, vers. 2.5.0 (Bouckaert et al., 2014). Each Markov 
chain Monte Carlo process was based on a run of 20 mil- 
lion generations, and genealogies were sampled every 1000 
generations. Of the sampled genealogies, 10% were dis- 
carded as burn-in. The divergence rate of 1.3% per million 
years was applied under the strict molecular clock model. 
The convergence of parameters was visually checked with 
effective sample sizes of estimates for all relevant param- 
eters over 200 with TRACER, vers. 1.7.1 (Rambaut et al., 
2018). The parameter sum (indicators.alltrees) that esti- 
mated the number of population change steps was used to 
check for the most likely numbers of demographic changes. 
Microsatellite analyses 
A PCR was performed to amplify each of the following 
7 SSR loci: LERI21, LERI26, LERI33, LERI34, LERI44, 
LERI50, and LERI63 (El Nagar et al., 2010). Polymerase 
chain reactions were carried out by using a PCR Thermal 
Cycler Dice Gradient, each with 10 pL of reaction solution 
containing 10 ng template DNA, 0.5 pM of each primer, 
0.2 mM dNTP mixture, 1x PCR buffer (10 mM Tris— 
HCl pH 8.3, 50 mM KCl, and 1.5 mM MgCl,), and 0.5 U 
