Ichthyol Res
DOI 10.1007/s10228-014-0408-9

FULL PAPER

Phylogeography of Blacktip Grouper, Epinephelus fasciatus(Perciformes: Serranidae), and influence of the Kuroshio Current
on cryptic lineages and genetic population structure
Kaoru Kuriiwa • Satoru N. Chiba •Hiroyuki Motomura • Keiichi Matsuura

Received: 4 December 2013 / Revised: 6 April 2014 / Accepted: 17 April 2014
Ó The Ichthyological Society of Japan 2014

Abstract To investigate the influence of the Kuroshio
Current on the high diversity of marine fishes in Japanesewaters, the intraspecific phylogeographic structure ofBlacktip Grouper, Epinephelus fasciatus, was determined.
The genetic analysis of E. fasciatus indicated three intraspecific mtDNA lineages representing different evolutionary histories: the first lineage differentiated in Japanesewaters during a long period of fluctuations of the ancient
Kuroshio Current, the second lineage, widely distributed in
the tropical western Pacific, was transported to Japanesewaters by the Kuroshio Current and the third lineage wasdistributed primarily around the Ogasawara (Bonin)
Islands. Present-day sympatric distributions of the threelineages, characterized by different ratios of such individuals at each geographic site, indicated a complex genetic
pattern that was classified into three demographic groups,
the dispersal and gene flows of which were strongly
influenced by the Kuroshio Current and factors such ascountercurrent and island arc. Genetic breaks in E. fasciatus populations were congruent with other fish faunalboundaries in Japanese waters.
Keywords Marine biogeography  Dispersal  Migration 
Pelagic larval duration

Introduction

Electronic supplementary material The online version of thisarticle (doi:10.1007/s10228-014-0408-9) contains supplementarymaterial, which is available to authorized users.
K. Kuriiwa  S. N. Chiba  K. Matsuura
National Museum of Nature and Science,
4-1-1 Amakubo, Tsukuba, Ibaraki 305-0005, Japan
S. N. Chibae-mail: chiba_s@kahaku.go.jp
K. Matsuurae-mail: matsuura@kahaku.go.jp
Present Address:
K. Kuriiwa (&)
Faculty of Science, Department of Chemistry, Biology
and Marine Science, University of the Ryukyus,
1 Senbaru, Nishihara, Okinawa 903-0213, Japane-mail: kaorukuriiwa@07.alumni.u-tokyo.ac.jp
H. Motomura
Kagoshima University Museum, 1-21-30 Korimoto,Kagoshima 890-0065, Japane-mail: motomura@kaum.kagoshima-u.ac.jp

The Kuroshio Current, one of the world’s major ocean
currents, flows from east of the Philippines to the Pacificcoast of southern mainland Japan, via Taiwan and thesouthern tip of the Ryukyu Islands. This current transports
tropical fishes from the Philippines and Taiwan to Japanesewaters resulting in overlapping distributions of tropical andwarm temperate fishes in the southern part of the latterregion (Nishimura 1992). Senou et al. (2006) indicated that
another important role of the current was to create aninvisible barrier to temperate fishes when they try to movesouthward. The Kuroshio Current and cold Oyashio Current (Kuril Current), extending southwards to northern
Japan, have resulted in a high diversity of marine organisms in Japanese waters (Fujikura et al. 2010). This
diversity has been the subject of frequent biogeographicalinvestigations (e.g., Jordan 1901; Tanaka 1931; Briggs
1974; Nishimura 1992; Nakabo 2002). Although several
faunal provinces in Japanese waters were defined on thebases of water temperature and the ocean currents, quantitative faunal data were not considered. In contrast, Senou

123

K. Kuriiwa et al.

Fig. 1 Color variations of Epinephelus fasciatus. a Kozu Island(collection locality), NSMT-P 96600 (registration number of voucher
specimen), 29.5 cm SL (standard length); b Shimokoshiki-jimaIsland, NSMT-P 95762, 31.4 cm SL; c Ogasawara Islands, NSMTP 95263, 24.6 cm SL; d Fukue Island, NSMT-P 91761, 29.0 cm SL;

e Tanegashima Island, KAUM-I 9069, 33.5 cm SL; f Yaeyama
Islands, NSMT-P 93854, 15.8 cm SL; g Kodakara-jima Island,
NSMT-P 91974, 18.1 cm SL; h Ogasawara Islands, NSMT-P 93541,
26.0 cm SL; i Ogasawara Islands, NSMT-P 95268, 23.2 cm SL

et al. (2006) and Matsuura and Senou (2012) conductedcluster analyses for the marine fish fauna of southern Japan,
which indicated two distinct faunal groups. In these faunalanalyses, however, the population sizes of fishes in eachgeographic site were not considered. In addition, althougha pair of closely related species distributed in each faunalgroup was considered, it is unclear whether cryptic genetic
breaks in widely distributed species around southern Japan
exist. A greater understanding of the fish diversity of thesefaunal groups requires detailed population or species levelbiogeographical analyses.
Phylogeography, a relatively recent biogeographicalapproach, is shedding light on the principles and processes governing the geographic distributions of genealogical lineages, especially those within and amongclosely related species (Avise 2000). The present study
focused on the intraspecific phylogeographic structure ofBlacktip Gouper, Epinephelus fasciatus, in Japanesewaters and the tropical western Pacific so as to determine how the Kuroshio Current may have influenced thisfish’s evolutionary history. Epinephelus fasciatus, one of
the most common and widespread species of the genus
in the tropical Indo-West Pacific, is also known fromsubtropical to temperate waters in Japan (Randall andHeemstra 1991). In addition, ecological characteristics ofthe species have been well documented in aquaculture
studies: (1) adult fishes spawn pelagic eggs, (2) pelagic
larval duration (PLD) is approximately two months, (3)larval settlement is at approximately 5 cm in total length(TL), (4) upon reaching 20 cm TL, most individuals are
reproductively mature, and (5) the species is protogynous, half of the subadult individuals changing sex tobecome male (Kawabe et al. 1997, 2000; Kawabe 2005;
Kawabe and Kohno 2009).
Our observations on E. fasciatus in Japanese waters
show differences in body coloration (Fig. 1) and maximum
body size based on geographic region. Most temperatewater individuals reach up to 50 cm TL, including thosearound the Izu Islands (populations 2–4 in Table 1 andFig. 2) and the main islands of Japan (5–13). These fishes
are characterized by bright red body coloration, whereas
those in tropical/subtropical waters around the Ogasawara(Bonin) Islands (1), the Ryukyu Islands (14–21) and further
south (22–27), have orange or whitish coloration and reachup to 35 cm TL. Many surveys using the phylogeographic
approach revealed that marine fishes generally show ashallower phylogeographic structure than that of freshwater/brackish-water fishes as a result of their continuous
habitat and relatively lengthy pelagic larval dispersal (Palumbi 1994; Grant and Bowen 1998; Avise 2000). However, the dual role of the Kuroshio Current and/or someother past events may cause discrete genetic structures in
E. fasciatus with phenotypic and quantitative variations
mentioned above. The phylogeographic structures andevolutionary history of E. fasciatus around Japan andadjacent waters are elucidated below.

123Phylogeography of Blacktip Grouper
Table 1 Collection localities, geographic population numbers andnumber of individuals examined of Epinephelus fasciatus
Locality

Geographic
population
numberNumber ofindividuals
examined

Japan, Ogasawara (Bonin) Islands

160

Izu Islands
Hachijo-jima Island

2

2

Kozu Island

3

32

Izu-oshima Island

4

31

Izu Peninsula

5

20

Mie Prefecture

6

20

Tokushima Prefecture

7

2

Kochi Prefecture

8

2Northern part of Miyazaki Pref.
(Kadogawa)

9

20Southern part of Miyazaki Pref.
(Nangou)

10

22

Fukue Island

11

20

Kagoshima Prefecture

121

Shimokoshiki-jima Island

13

16

Tanegashima Island

14

19

Yaku-shima Island

15

22

Kuchino-shima Island

16

10

Gaja-jima Island

17

20

Kodakara-jima Island
Amami-oshima Island

18
19

16
10

Okinawa Island

20

25

Yaeyama Islands

2130

Taiwan, East coast (Su-ao, Hualien,
Cheng-gon, and Fugang)

22

25Southern coast (Hengchun)

23

20

West coast (Wuchi and Penghu)

24

12Vietnam, Nha Trang

25

35

Malaysia, Kota Kinabalu

26

26

Indonesia, Bitung (Northern tip of the
Sulawesi Island)

27

23

Main islands of Japan

Ryukyu Islands

Materials and methods
Fish samples. In total, 541 individuals of Epinephelus
fasciatus were used in this study, comprising specimensfrom Japan (400), Taiwan (57), Vietnam (35), Malaysia(26) and Indonesia (23) (Fig. 2; Table 1). Thirteen additional species of Epinephelus were selected as outgroupspecies [Electric Supplementary Material (ESM) Table
S1], based on phylogenetic analyses and a revised classification of the tribe Epinephelini (Craig and Hastings2007). Species identification and nomenclature followed
Randall and Heemstra (1991) and Senou (2002).

Specimens were obtained by spear-fishing, line-fishing orfrom fish markets between 2007–2011. Voucher specimens
were deposited in the National Museum of Nature andScience (NSMT - formerly National Science Museum,
Tokyo), the Kagoshima University Museum (KAUM), theFisheries Research Laboratory, Mie University (FRLM)
and Division of Fisheries Sciences, Faculty of Agriculture,
Miyazaki University (MUFS). All tissue samples andextracted DNA were deposited in NSMT. Specimen details
and DNA data for the fishes used in this study are given in
Table S1. The map with submarine topographies was created by the Generic Mapping Tools (GMT) ver. 4.5.7(University of Hawaii: http://gmt.soest.hawaii.edu/).
DNA extraction, PCR amplification and sequencing.
Total genomic DNA was extracted from the tissues using
the GentraÒ PuregeneTM Tissue Kit (QIAGEN), following
the manufacturer’s protocol. A segment of approximately
2.3 kbp of mitochondrial DNA (mtDNA) was amplified by
the polymerase chain reaction (PCR) with a pair of primers. The primer sequences were Epi-cyF: 5’-ATG GCH
AGC CTH CGH AAA ACN CA-3’ as the forward primer
and Simt12-Rb: 5’-TGC RGA TAC TTG CAT GTG TAA3’ as the reverse primer. The former was designed withreference to the cytochrome b (cyt b) gene sequences of E.
fasciatus (DDBJ/EMBL/GenBank accession number
AY786423), Epinephelus merra (AY786424), Epinephelus
awoara (AB108494), Epinephelus tauvina (EF456003),
Epinephelus septemfasciatus (AY786425), Cephalopholis
urodeta (AY786426) and Variola louti (AY786428). The
latter was the same as Simt12-Fb, described in Kuriiwa
et al. (2007). The amplified region using this primer pair
was from the 5’ end of the cyt b gene to upstream of the 5’of the 12S rRNA gene, thereby including completesequences of the cyt b gene, tRNA-Threonine (tRNA-Thr),tRNA-Proline (tRNA-Pro) and control region (CR).
PCR was carried out with a 10 ll reaction volumecontaining 1.0 ll of Ex Taq Buffer, 0.8 ll dNTP mixture,
0.1 ll of Ex TaqÒ DNA polymerase (TaKaRa), 7.2 llsterile distilled water, 0.2 ll of each primer (10 lM) and
0.5 ll of template DNA. Amplification parameters were 35
cycles of denaturation at 94 °C for 30 s, annealing at 52 °C
for 30 s and extension at 72 °C for 60 s. Double-stranded
PCR products, purified using the ExoSAP-IT (USB), were
subsequently used for direct cycle sequencing with the
BigDyeÒ Terminator Cycle Sequencing Kit ver. 3.1(Applied Biosystems), following the manufacturer’s protocol. Six primers were used for sequencing the PCR
products: those described above for PCR and four newly
designed (Epi-cyF2: 5’-TAT TCC TAA TAC CAG CAG
C-3’ and Epi-CRF: 5’-AAG ACA AGG ACT YCT TGA
AGG TC-3’ as forward primers, and Epi-proR: 5’-GTTTAA TTT AGA ATT CTA GCT TTG G-3’ and Epi-cyR:5’-TTG TCG GCA TCA GAR TTA AGG CC-3’ as reverse

123

K. Kuriiwa et al.

Fig. 2 Collection localities of Epinephelus fasciatus. For details of geographic population numbers and number of individuals examined, see
Table 1primers). Dye-labeled fragments were analyzed on a model
3130xl Genetic Analyzer (Applied Biosystems). The DNA
sequences were edited and aligned with ATGC ver. 5.0.1(GENETYX) and MEGA ver. 5.1 (Tamura et al. 2011). All
sequences are available from DDBJ/EMBL/GenBank
accession numbers AB705627–AB706181 and AB829908–
AB829938 (ESM Table S1).
Phylogenetic and haplotype network analyses. A maximum likelihood (ML) tree was inferred by LIKELIHOODRATCHET (Vos 2003), using all haplotypes of the wholesegment sequences of the cyt b, tRNA-Thr, tRNA-Pro andCR for E. fasciatus. Owing to the long length of the CR
sequences for 10 of the 13 outgroups, we succeeded indetermining the sequences of only the other three outgroups being adequate for all analyses below. Haplotypes
were detected by the program TCS (see outline in haplotype network analysis). Initially, 1,000 trees were generated by the program TNT ver. 1.1 (Goloboff et al. 2008)
with 25 % random site-upweighted data, created by theparallel execution program PHYLOGEARS ver. 1.5(Tanabe 2008). Subsequently, these trees were subjected toa heuristic search using TREEFINDER ver. Oct. 2008
(Jobb et al. 2004; Jobb 2008), with application of the TN93
(Tamura and Nei 1993) ? I ? C model, which was selected
as the best-fit model of nucleotide substitution using KAKUSAN ver. 4.0 (Tanabe 2007), on the basis of the partition- and codon-proportional model with Akaike
Information Criterion (AIC, Akaike 1973). To account forthe difference rate among genes and the codon positions,123

we divided the whole segment sequences into the following
six partitions; the 1st, 2nd, and 3rd codon sites of the cyt b,
tRNA-Thr, tRNA-Pro and CR. Support for the internalbranches of the eventuating tree was assessed using bootstrap resampling with 1,000 replicates by TREEFINDER
with PHYLOGEARS.
Phylogenetic trees were also reconstructed by Bayesian
phylogenetic analysis, conducted with MRBAYES ver.
3.2.1 (Ronquist et al. 2012), using the same dataset, partitions and AIC parameters. This program can use different nucleotide substitution model on each partition, thus
we set the HKY85 (Hasegawa et al. 1985) ? I, GTR(Tavaré 1986) ? C, K80 (Kimura 1980) ? I,
K80_homogeneous, HKY85 ? C and GTR ? C model,
which was selected as the best-fit model of nucleotide
substitution using KAKUSAN, on the first, second, andthird codon sites of the cyt b, tRNA-Thr, tRNA-Pro andCR, respectively. Two independent runs of 12,000,000
generations were conducted in each analysis, with tree
sampling every 100 generations and burn-in after 12,001
trees. The remaining 108,000 trees from the two independent analyses were used to generate a 50 % majority
rule consensus tree, with the percentage of samples
recovering any particular clade representing the posterior
probability.
On the basis of the average nucleotide substitution rateof the cyt b gene and CR, 2 % and 10 % per million years,
respectively, between lineages of bony fishes (Bowen et al.
2006b), the divergence time among genetic clades wasPhylogeography of Blacktip Grouper
Table 2 Genetic variations, results of the neutrality test and demographic indices for sequences of mtDNA cytochrome b gene (1,053 bp),
control region (872 bp) and whole segment (2,068 bp) within each intraspecific lineage of Epinephelus fasciatus
Lineage

n

Cyt b

CR

Whole segment

H

H

H

UST

p 9 100h

SSD (P)

Ragged (P)

A

129

19

44

57

0.873

0.931±0.015

0.324±0.170

0.014 (0.408)

0.015 (0.574)

B

386

139

224

291

0.873

0.995±0.001

0.440±0.224

0.002 (0.710)

0.003 (0.890)

C

26

13

21

23

0.872

0.991±0.013

0.358±0.193

0.004 (0.547)

0.013 (0.648)

Total

541

171

289

371

0.872

0.994±0.001

1.594±0.770--(Continued)
Lineage

Whole segment (continued)
Tajima’s Ds

Fu’s FS

h0

h1

h1/h0

A

-1.444** (0.040)-24.914* (0.000)

11.0120.002

12.4026201.000B-2.415* (0.000)-24.038* (0.002)

10.6190.009

36.8024089.111

C-2.076* (0.000)-14.144* (0.000)

1.002

78.750

Total

-

-

-

-

5.439
-

78.593-

Cyt b cytochrome b; CR control region; Whole whole segment of the cyt b, tRNA-Threonine, tRNA-Proline and CR; n number of individuals
examined; H number of haplotypes; h haplotype diversity; p nucleotide diversity; SSD sum of squared deviations; Ragged raggedness index; sunits of mutational time; h0 population size before expansion; h1 population size after expansion
Asterisks (*) and (**) indicate significance at [ 99 % and [95 % confidence interval, respectively

Table 3 The average p-distance (%) of the cytochrome b gene (1,053 bp), control region (872 bp) and whole segment (2,068 bp) sequences
within species and each lineage, and between lineages of Epinephelus fasciatus
Within species (range)

Within lineage (range)

Between lineages (range)

A

B

C

A and B

A and C

B and C

Cyt b
CR

1.2 (0–4.7)
2.0 (0–8.0)

0.2 (0–0.7)
0.6 (0–1.7)

0.4 (0–1.1)
0.6 (0–1.9)

0.2 (0–0.8)
0.6 (0–2.7)

1.9 (1.5–2.5)
3.1 (2.3–4.5)

3.8 (3.5–4.2)
7.3 (6.6–8.0)

3.8 (3.4–4.7)
6.8 (5.9–8.0)

Whole

1.6 (0–5.9)

0.3 (0–0.8)

0.4 (0–1.2)

0.3 (0–1.3)

2.6 (2.2–3.2)

5.6 (5.2–5.9)

4.9 (4.4–5.5)

Cyt b cytochrome b; CR control region; Whole whole segment of the cyt b, tRNA-Threonine, tRNA-Proline and CR

roughly estimated. Subsequently, the molecular clock
hypothesis was tested with the maximum likelihood
method (likelihood ratio test; LRT) using the whole segment sequences for E. fasciatus by the BASEML programof PAML ver. 4.7 (Yang 2007). To avoid the computational burden, 50 haplotypes were randomly selected from371 haplotypes. These consist of Lineages A, B and C. The
LRT was carried out on the basis of these data sets andrepeated three times independently.
To clarify the historical process of mtDNA haplotypes, astatistical parsimony network in each clade of the phylogenetic tree was constructed by TCS ver. 1.21 (Clement
et al. 2000), which employs the method of Templeton et al.
(1992). The networks of both the CR and whole segment
sequences were too reticulate to detect any structures (datanot shown), so that we used the cyt b gene sequences for all
individuals of E. fasciatus. The networks were connectedusing a maximum parsimony algorithm, with a statisticalconfidence of 95 %. The ancestral haplotype was also

identified using TCS by the method of Castelloe andTempleton (1994), which estimates outgroup-weights
based on haplotype frequency and connectivity, the haplotype with the highest outgroup-weight most likely being
the oldest.
Demographic analyses and population structure. The
program ARLEQUIN ver. 3.11 (Excoffier et al. 2005) wasused to generate haplotype diversity (h) and nucleotide
diversity (p) of the whole segment sequences within
genetic clades. This program was also used for the following analyses and tests. Theoretical distributions undersimulated constant population size and a sudden expansion
model (Rogers and Harpending 1992) were compared to
the goodness-of-fit of the observed data with the sum ofsquared deviations (SSD) and raggedness indices (Herpending 1994). Neutrality of the sequence variation of eachgenetic clade was verified using Tajima’s D (Tajima 1989)
and Fu’s FS (Fu 1997) tests. To estimate demographic
parameters, pairwise mismatch distributions (Rogers and123

K. Kuriiwa et al.

Fig. 3 Phylogenetic tree and geographic distribution of intraspecificlineages for Epinephelus fasciatus. a Maximum likelihood (ML) tree
based on whole segment sequences of mtDNA cytochrome b gene,tRNA-Threonine, tRNA-Proline and control region from 371 haplotypes using the TN93? I ? C model of sequence evolution on thebasis of the partition- and codon-proportional model with Akaike
Information Criterion. Numerals above nodes indicate bootstrap
values in the ML (left) and posterior probabilities in the Bayesiantrees (right) from 1,000 replicates. Solid and dotted line scale
indicates 0.01 % (E. fasciatus) and 0.1 % (outgroups) substitution persite, respectively; b geographic distribution of intraspecific lineagesof E. fasciatus. Each color code of circle refers to three lineages in thephylogenetic tree (a). Circle size proportional to number of individuals. For details of geographic population numbers beside open
circles and number of individuals examined, see Table 1 and ESM
Table S2

Harpending 1992) were calculated within genetic clades.
The times and magnitudes of inferred sudden expansions
were determined by calculating s (assuming s = 2 ut,
where u is the mutation rate per generation for the entire
region of DNA under study and t is the time in generations), h0 and h1 (assuming h0 = 2 N0l and h1 = 2 N1l,
where N0 and N1 are the population sizes before and after
expansion, and l is the mutation rate per nucleotide), using
a set of 1,000 parametric bootstrap replicates (Rogers andHarpending 1992; Rogers 1995; Schneider and Excoffier
1999; Excoffier 2004).
The Bayesian skyline plot (BSP; Drummond et al. 2005)
for each genetic clade was carried out to clarify past
dynamics in the female effective population size (Nef)through time, using BEAST ver. 1.7.5 (Drummond andRambaut 2007). The best-fit model of nucleotide substitution for whole segment sequences was selected using
KAKUSAN. Each Markov chain Monte Carlo (MCMC)
process was based on a run of 10,000,000 generations, andgenealogies were sampled every 1,000 generations. Tenpercent of sampled genealogies were discarded as burn-in.
On the basis of the harmonic mean of nucleotide divergence rates of the cyt b gene and CR between lineages ofbony fishes (Bowen et al. 2006b), nucleotide substitution
rate of 3.3 % per million years was applied under the strict
molecular clock model based on the results of the LRTanalyzed by BASEML. Convergence for the posteriormeans of all parameters was assessed using TRACER ver.
1.5 (Rambaut and Drummond 2009) by checking
Nef [ 200.
The genetic structure for each of 23 geographic populations, which included more than 10 individuals, wasassessed by an analysis of molecular variance AMOVA(Excoffier et al. 1992), as implemented in ARLEQUIN.
The diversity measure U-statistics values were calculatedusing the whole segment sequences. The criterion for
significant differences of pairwise UST (genetic differences between populations) values from zero among the23 geographic populations was tested by a nonparametricpermutation test (Rice 1989). Population grouping based
upon genetic variability was performed by SAMOVAver. 1.0 (Dupanloup et al. 2002), which pooled populations into alternative structures by maximizing UCT(genetic differences among the population groups)
values.

123

Results
Genetic variation. DNA segments of 2,061–2,068 bp weresuccessfully sequenced for all 541 individuals of Epinephelus fasciatus. These included partial sequences of thecyt b gene (1,053 bp) and complete sequences of thetRNA-Thr (73 bp), tRNA-Pro (69–70 bp) and CR

Phylogeography of Blacktip Grouper

Fig. 4 Haplotype networks of the cytochrome b gene sequences fromall 541 individuals of intraspecific lineages for Epinephelus fasciatus.
Circles connected to one another represent mtDNA haplotypes ofcytochrome b gene sequences (1,053 bp). Circle size is proportional

to haplotype frequency. Small open circles represent putative(missing) haplotypes not detected in this study. Color codes of eachcircle refer to geographic populations in map (above). Solid arrow
indicates highest outgroup-weight haplotype in each lineage(866–872 bp). A total of 170 (cyt b), 8 (tRNA-Thr), 6(tRNA-Pro), 229 (CR) and 413 (whole segment) nucleotide
positions varied within the species, such variations defining
171 (cyt b), 9 (tRNA-Thr), 6 (tRNA-Pro), 289 (CR) and
371 (whole segment) haplotypes (Table 2). The same
regions of DNA segments for three outgroups (E. retouti,
E. akaara and E. awoara) were also determined(2,184–2,293 bp), only the cyt b gene sequences being
determined for the remaining ten due to the long length oftheir CR sequences (ESM Table S1). The average p-distance (uncorrected sequence divergence) of the cyt b, CR
and whole segment sequences within E. fasciatus areshown in Table 3. That of the cyt b was 14.9 % (range
9.6–16.3 %) between E. fasciatus and 13 outgroups, and123

K. Kuriiwa et al.
Fig. 5 Mismatch distributions
and Bayesian skyline plots(BSPs) from all 541 individuals
of intraspecific lineages for
Epinephelus fasciatus.
a Observed (vertical bars) and
expected (solid lines) mismatch
distribution of pairwise
sequence divergence for wholesegment sequences of mtDNA
cytochrome b gene, tRNAThreonine, tRNA-Proline andcontrol region. Horizontal andvertical axes indicate genetic
differences and frequency,
respectively. Each color coderefers to each lineage in thephylogenetic tree (Fig. 3a);b BSPs from the same dataset ofmismatch distributions.
Horizontal axis represents timein million years ago and verticalaxis representing the female
effective population size (Nef) 9
generation interval (10-6
million years ago) is given on alogarithmic scale. Solid anddotted lines represent mean
Nef 9 10-6 and 95 % highest
probability density (HPD)
intervals, respectively. Inserted
gray regions indicate the Last
Glacial Maximum (LGM; ca.
0.01–0.02 million years ago)

14.1 % (6.0–16.5 %) among the 13 outgroups. That of the
CR and whole segment sequences between E. fasciatus andthree outgroups was 20.3 % (17.2–23.6 %) and 16.5 %(14.4–18.4 %), and among three outgroups was 23.9 %(20.9–25.5 %) and 17.0 % (13.5–19.0 %), respectively.
Phylogenetic relationships. The phylogenetic analyses
using the 371 haplotypes of the whole segment sequences
revealed that E. fasciatus was divided into three highly
supported genetic clades (Fig. 3a; ESM Fig. S1), referredto here as intraspecific lineages A, B and C. Both the ML
and Bayesian analyses produced an identical topology.
Each lineage consisted of 57 haplotypes for 129 individuals(Lineage A), 291 for 386 (B) and 23 for 26 (C), from many
geographic populations (Fig. 3b; Table 2; ESM Table S2).
The average p-distance of the cyt b, CR and whole segment
sequences within a single lineage and between lineages arelisted in Table 3. The divergence time of lineages A/B and
AB/C of E. fasciatus was roughly estimated as having been

123

0.95 and 1.90 (cyt b), and 0.31 and 0.71 (CR) million years
ago (MYA), respectively. The strict molecular clock
hypothesis could not be rejected in every trial of the LRTanalyses (P [ 0.05; ESM Table S3).
Haplotype network. Haplotype networks of the cytb gene sequences for the three lineages are shown in
Fig. 4. Those for Lineage A exhibited two dominant
haplotypes with several rare haplotypes that differed byless than three mutational steps. The two dominant haplotypes were connected by two mutational steps. The
network for Lineage B was more complex, consisting offour dominant haplotypes, connected by one to two
mutational steps, and numerous rare haplotypes that differed by less than six mutational steps. Although thenumber of individuals was insufficient for determining adefinite structure, the Lineage C network included a single
dominant haplotype with several rare (one highly differentiated) haplotypes.

Phylogeography of Blacktip Grouper
Fig. 6 Phylogeographic
provinces of Epinephelus
fasciatus. a Frequency of
intraspecific lineages in 23
geographic populations (above)and phylogeographic provinces
indicated by populationgrouping (below). Geographic
population numbers (oblique)and color codes refer to Table 1
and Fig. 3, respectively. K value
indicates number of groups for
maximizing genetic differencesamong population groups, UCT
values. Roman numerals
indicate population group
numbers. For details of
population grouping, see ESM
Table S3; b bold and thinarrows represent routes of theKuroshio Current and Kuroshio
Countercurrent, respectively.
Dotted lines indicate
phylogeographic boundaries
with reference tophylogeographic provinces of E.
fasciatus (a). Geographic
population numbers inside open
circles refer to Table 1

Demographic analyses. The demographic analyses
using the whole segment sequences showed a historicaldemography of each lineage. Although the haplotype
diversity (h) was very high (h = 0.931–0.995), nucleotide
diversity (p) was low (p = 0.003–0.004) in all lineages(Table 2), suggesting sudden expansions after bottleneck
effects. Neutrality tests also showed evidence of sudden
expansions in all lineages; insignificant SSD values andraggedness indices (P [ 0.05), and significant Tajima’s
D and Fu’s FS values (P \ 0.05) occurred for all lineages(Table 2).
Although mismatch distributions in Lineage B were
unimodal and consistent with an expanding population
size, Lineages A and C represented bimodal distributions(Fig. 5a). In the bimodal distribution of Lineage A, twopeaks were connected with a continual degree of

pairwise differences. In that of Lineage C, the lower
peak was caused by a highly differentiated individual,
reflecting the phylogenetic tree and haplotype network(Figs. 3a, 4, 5a).
The BSPs showed that the population expansion
dynamics in each lineage at approximately 0.05 (LineageA), 0.085 (B) and 0.115 (C) MYA of the coalescent time
for the most recent common ancestor (MRCA), respectively, from the analysis with the strict molecular clock
model (Fig. 5b). The BSPs also showed that the Nef clearly
departed from the constant population size model in eachlineage (Fig. 5b). In Lineage A, although the BSP did notshow departure from the constant size model before theLast Glacial Maximum (LGM; ca. 0.01–0.02 MYA), theNef expanded notably after the LGM. The BSP in Lineage
B showed gradual increase before the LGM but notable123

K. Kuriiwa et al.

Fig. 7 Color-coded visible columns of pairwise UST values among
23 geographic populations of Epinephelus fasciatus. Each color coderefers to degree of UST values (right). Asterisks indicate significanceat 95 % confidence interval. Geographic population numbers(oblique) refer to Table 1. Dotted lines indicate phylogeographic
boundaries of E. fasciatus (Fig. 6a). Raw data of pairwise UST and
P values refer to ESM Table S5. The style of this figure was referredto that of fig. 2 given by Hickey et al. (2009)

after the LGM. On the contrary, that in Lineage C wasconsistently and gradually increasing, being not related to
the LGM.
Demographic indices using the whole segment sequences indicated the time and magnitude of expansion of eachlineage (Table 2). The observed s value of Lineage C was
smaller than in Lineages A and B, indicating that the timeat which expansion started was later in Lineage C than in
the other two lineages. Although the estimated effective
population size before expansion (h0) in Lineage C wasmuch greater than zero, whereas those in Lineages A and B
equaled approximately zero, the paces of expansion (h1/h0)
in Lineages A and B were greater than in Lineage C. A
comparison of the h1 values suggested that the present
population size of Lineage C is the largest, that of Lineage
A being the smallest.
Population structure. Each geographic populationconsisted of individuals of the three genetic lineages, albeit
in different ratios (Figs. 3b, 6a; ESM Table S2). Twentythree geographic populations were divided into threedemographic groups (Groups I, II and III) by SAMOVA
analyses (Fig. 6a; ESM Table S3), in which the K value(number of groups) = 3 was selected by maximizing UCT
values (UCT = 0.342, USC = 0.003, UST = 0.344). Group I
consisted of a population from the Ogasawara Islands(population-1), including individuals of all three lineages,
each in an approximately equal ratio. Group II included bothof the geographic populations from the Izu Islands (3 and 4)
and most from the main islands of Japan (5, 6, 9, 11 and 13),
including primarily individuals of Lineage A, but also thoseof Lineage B. Group III consisted of a population from themain islands of Japan (10) (dominated by individuals of123

Lineage B but also including those of Lineage A), andpopulations from the Ryukyu Islands (14–21) and all countries sampled south of Japan (22–27), all of which comprised
individuals of Lineage B only. Mismatch distributions of 23
geographic populations and three demographic groups areshown in ESM Figs. S2 and S3, respectively. In those, thenumber of peaks for pairwise differences depended on thenumber of genetic lineages being included in each geographic population and demographic groups.
The three demographic groups were also supported by the
results of the population pairwise UST analyses (Fig. 7; ESM
Table S4). Most of the significant differentiation foundoccurred among groups (UST = 0.071–0.558; P \ 0.05),
rather than within groups, the exception being population10. Although this population was included in Group III, asignificant differentiation was found not only between populations of Group I (UST = 0.117; P \ 0.01) and II(UST = 0.107–0.218; P \ 0.05), but also between otherpopulations of Group III (UST = 0.067–0.159; P \ 0.05).

Discussion
Evolutionary history of E. fasciatus. The molecularphylogenetic and demographic analyses indicated that
Epinephelus fasciatus distributed in Japanese waters andthe tropical western Pacific is genetically clustered into
three cryptic intraspecific lineages (Fig. 3a: Lineages A, B
and C), which exhibit different evolutionary histories.
Lineage A (termed the ‘‘Japanese lineage’’) probably
became differentiated in Japanese waters in the early
middle Pleistocene (ca. 0.13–2.59 MYA) under the influence of long-period fluctuations in the route of the Kuroshio Current (Fig. 3a). During the Pleistocene (ca.
0.01–2.59 MYA), land bridges connecting the main islandsof Japan to continental Asia, and the central-southern Ryukyu Islands to continental Asia via Taiwan, appeared anddisappeared several times (Ota 1998; Kimura 2000). When
land bridges appeared, especially during the LGM (ca.
0.01–0.02 MYA), the ancient Kuroshio Current turned tothe east at the southern end of the Ryukyu Islands (Ujiié
and Ujiié 1999; Ujiié et al. 2003). When land bridges
disappeared, however, the ancient Kuroshio Current may
have flowed northward via the southern tip of the Ryukyu
Islands and Taiwan, a route similar to that of the presentday. The continual bimodal distribution in a mismatchanalysis and the population expansion after the LGM in the
BSP of Lineage A (Fig. 5a, b) suggested that some populations dispersed and were isolated when land bridges
appeared, and subsequently experienced the secondary
contact when they disappeared after the LGM.
Although individuals of Lineage B, here termed the‘‘northwestern Pacific lineage’’, were found in allPhylogeography of Blacktip Grouper

geographic populations and those of Lineages A and B
were sympatric in Japanese waters (Figs. 3b, 6a), thenorthern individuals of Lineage B were primarily found in
the Ryukyu Islands, their distribution extending northward
to the coast of the main islands of Japan through transportation by the Kuroshio Current. An alternative scenario,
based on differences in adaptive abilities, such as tolerance
to cold temperatures, might also be considered. Individualsof Lineages A and B were in greater proportion in thetemperate region (main islands of Japan) with the opposite
being true for the tropical/subtropical region (Ryukyu
Islands). A comparison of the ratios of these two lineages
between populations 9 vs 10, 11 vs 13, 9 vs 11 and 10 vs
13, showed that of Lineage B to be higher than that ofLineage A in the area strongly influenced by the Kuroshio
Current and thereby characterized by a higher temperature(Fig. 6). However, in the subtropical Ogasawara Islands(population-1), individuals of these two lineages were
found in almost equal ratio. Given these data, the northward dispersal of Lineage B by the Kuroshio Current is areasonable hypothesis. The BSP of Lineage B suggested
that the ancient Kuroshio Current was mainly influenced
only in the Ryukyu Islands and the northernmost individuals of Lineage B were distributed in the Ryukyu Islands
until the LGM. However, the route of the Kuroshio Current
changed to the present flow after the LGM, expanding the
distribution and the population size of E. fasciatus(Fig. 5b). The distributions of Lineages A and B extended
along the route of the Kuroshio Current (Fig. 3b), and all
dominant haplotypes in the haplotype networks for these
two lineages occur in individuals from most geographic
populations (Fig. 4). This implies that most of the dominant haplotypes existed before expansion of the populations with large geographic-scale gene flow occurring after
expansion due to current-assisted transport of individuals.
The evolutionary history of Lineage C most likely differs from those of the other two lineages due to the geographic distribution of the former. The Ogasawara Islands,
where most of the individuals of Lineage C were found, aretypical oceanic islands and not strongly influenced by the
Kuroshio Current. The BSP of Lineage C suggested that
not only the present day current but also the ancient one did
not influence the islands (Fig. 5b). In addition, the center of
distribution of Lineage C appears to be outside the areafrom which the present specimens were sampled, althoughmany individuals of this lineage were found in the Ogasawara Islands (Fig. 3b; ESM Table S2). This reasoning
was based on the large value for the estimated population
size of Lineage C after expansion, in spite of the lineage
being represented by the lowest number of individuals
sampled (Table 2). The most differentiated individual(Figs. 3a, 4, 5a) may also support the hypothesis. We
suspect that individuals of Lineage C are mainly distributedin the Mariana Islands, located on the eastern edge of thePhilippine Sea Plate. Not only did Springer (1982) proposea pattern of distribution and endemicity of shorefish species
associated with that tectonic plate, but also Randall et al.
(1997) demonstrated an affinity of fish faunas between the
Ogasawara and Mariana islands. Specimens from the Mariana Islands are required for clarification of the distribution pattern of Lineage C.
Interestingly, the three genetic lineages are not distinguishable by body coloration or maximum body size.
Individuals a–c, d–f and g–i in Fig. 1, belonging, respectively, to Lineages A, B and C, clearly illustrated that eachgenetic lineage includes all color variations. We note thatorange or whitish individuals were found in tropical/subtropical regions and reddish ones in temperate regions
suggesting that coloration may be related to habitat factors
rather than mitochondrial DNA lineages. Aquacultured
individuals of E. fasciatus in the Ogasawara Islands provided a clue with hatchery-reared juveniles being darkbrown (ESM Fig. S2), changing to red-, orange- or whitetinged within a year of their release (unpublished observation by the Tokyo Metropolitan Ogasawara Fisheries
Research Center). This suggests that coloration of E. fasciatus may be influenced by diet.
Maximum body size (ESM Table S6) may also be
related to habitat factors, most likely water temperature andthe length and/or times of reproduction. Recent studies ofhatchery-reared E. fasciatus revealed that the length of thespawning period in a subtropical region (Ogasawara
Islands) was more than twice that in temperate waters(main islands of Japan) (Kawabe et al. 2000). In addition,
the spawning frequency of individuals in the Ogasawara
Islands is possibly more than once per season (unpublishedobservation by the Tokyo Metropolitan Ogasawara Fisheries Research Center). Moreover, the average diameter ofpelagic eggs and the total length of larval E. fasciatus havea significant negative correlation with increasing watertemperature (Kawabe et al. 2000; Kawabe 2005). Individuals living in temperate regions may, therefore, grow larger
than those in tropical/subtropical regions.
The sympatric occurrences of the three DNA lineages atthe collection sites suggested that individuals can successfully interbreed irrespective of their genetic lineages.
An analysis of 32 hatchery-reared E. fasciatus individuals
from the Fisheries Cooperative Association of Chichi-jimaIsland (Ogasawara Islands), which were randomly selectedfrom clutches spawned in 2004–2007, found all three lineages, in the ratio A:B:C = 13:18:1. In that fisheries
cooperative, spawned eggs are obtained from several
hundred, multi-generational adult fish, implying an absenceof genetic incompatibility among the DNA lineages.
Genetic population structure of E. fasciatus. The
present study of E. fasciatus from many localities in123

K. Kuriiwa et al.

Japanese waters and the tropical western Pacific revealed
the distinctive genetic population structure of the species,
its uniqueness lying in the occurrence of individuals of
three mitochondrial DNA lineages having sympatric distributions with different lineage ratios in each geographic
population, hence three demographic groups (I, II and III)(Figs. 3b, 6a). Such a complex structure in a species ofpelagic larval duration (PLD) approximately two months
comes as a surprise because the larvae may be able todisperse over a long distance. Many previous studies haveinvestigated the relationships between PLD and genetic
population structure in marine fishes, those having a shortor no PLD generally exhibiting an obvious genetic structure (e.g., Riginos and Nachman 2001), whereas those witha lengthy PLD may have a non-significant genetic structure. The population structure in marine fishes may be
determined not only by PLD but also larval retention andself-recruitment (Jones et al. 1999; Swearer et al. 1999;
Almany et al. 2007), adult habitat preference (Rocha et al.
2002), historical vicariance (Bernardi et al. 2003; Shirai
et al. 2006; Liu et al. 2007), dispersal ability and endemic
distribution (Bowen et al. 2006a; Eble et al. 2009), ocean
currents (Rocha et al. 2008) and genetic isolation-by-distance (Hickey et al. 2009). In the case of E. fasciatus in
Japanese waters, the dual roles of the Kuroshio Current,
those of a transport system and an unseen marine barrier,
may contribute concurrently in determining the population
structure. The phylogeographic pattern suggested that
individuals of Lineage B can disperse from the Ryukyu
Islands and tropical western Pacific (Group III) to the mainislands of Japan (II), transported by the Kuroshio Current.
However, those of Lineage A cannot be transported from
north to south, due to the Kuroshio Current barrier (Fig. 6a,b).
It has been well established that the Kuroshio Current
fluctuates greatly both north and south of the central island(Honshu) of Japan being classified into five types (Fig. 6b;data from the Hydrographic and Oceanographic Department, Japan Coast Guard). Such fluctuations provide forboth dispersal and migration of marine fishes along thecoast of the central island of Japan and the Izu Islands, andsometimes also the Ogasawara Islands (Senou et al. 2006).
Although the distribution of individuals of Lineages A and
B in Fig. 6a strongly supports a hypothesis of fish dispersal
from north to south (from Group II to I), dispersal in the
opposite direction (from I to II) was not indicated by the
distribution of Lineage C. Several endemic species in the
Ogasawara Islands including Chaetodon daedalma and
Scarus obishime have been often found in the southern partof the Izu Islands (Senou et al. 2002, 2006), supporting an
hypothesis of fish dispersal from south to north. However,opportunities for northward and southward dispersal differ:
whereas the Kuroshio Current may play a positive role in123

southward transport along the Izu-Ogasawara (Bonin) Arc(see below), it may act as a weak barrier for northerly
dispersal. This hypothesis cannot be further investigatedhere due to an insufficient number of individuals obtained
from Hachijo-jima Island (population-2) (southern Izu
Islands).
In addition to the dual roles of the Kuroshio Current
mentioned above, a geographic factor is also thought tocontribute to the population structure between the Izu andOgasawara islands, which together constitute the IzuOgasawara Arc, expanding for approximately 1,200 km
from Honshu Island to the Volcano Islands (Fig. 6b).
Generally, dispersed larvae cannot survive unless theyfind shallow coastal or island waters to settle, whichsuggests the ‘‘stepping stone theory’’ (see ‘‘stepping stone
model’’ of population structure in Kimura and Weiss
1964). It is, therefore, reasonable to consider that fishlarvae transported by the Kuroshio Current can survive inthis island arc. Moreover, not only larvae but also adult
fishes may be able to migrate along this island arc.
Several exceptional individuals representing Lineages A
and C in the Ryukyu Islands population (Figs. 3b, 6a)
indicated that larval dispersal of E. fasciatus rarely occurs
from the Ogasawara Islands to the Ryukyu Islands. Despitea weak westerly current (Kuroshio Countercurrent), its
unsettled nature and lack of major island ‘‘steppingstones’’, except for the Daito (Borodino) Islands comprising several very small islands in the extensive area covering some 1,400 km between the two island groups(Fig. 6b), make for difficult westward movement for mostfishes. However, a few species of coral reef fishes havebeen reported as having been transported from the Ogasawara Islands to the Daito and/or Ryukyu islands: Chaetodon daedalma (see Yoshigou 2004), Hemitaurichthys
thompsoni (see Shimada 2002) and Ctenochaetus hawaiiensis (see Senou et al. 2003). We, therefore, consider thatE. fasciatus, having a lengthy PLD, may disperse accidentally by means of the Kuroshio Countercurrent.
Marine biogeographic boundaries with reference tofishes in Japanese waters. An additional finding of thisstudy was that the genetic population breaks in E. fasciatus(Fig. 6a, b) coincided with the marine biogeographic
boundaries (with reference to the fish fauna in Japanesewaters), proposed by Matsuura and Senou (2012). Cluster
analyses of fish fauna by them revealed that 12 regions insouthern Japan were grouped into two clusters: the firstcomprising islands located in the Ryukyu Islands, and theother comprising the main islands of Japan, and the Izu andOgasawara islands. In the latter cluster, a weak faunal
break was found between the Izu and Ogasawara islands.
Acknowledgements We express our sincere gratitude to the following persons and organizations for helping in the collection ofPhylogeography of Blacktip Grouper
fishes: Y. Iwatsuki (MUFS); K.T. Shao (Research Center of Biodiversity, Academia Sinica, Taiwan); H.C. Ho (National Museum ofMarine Biology and Aquarium, Taiwan); W.C. Jiang and colleagues(Eastern Marine Biology Research Center, Fisheries Research Institute, Taiwan); S. Kimura (FRLM); G. Ogihara and M. Matsunuma(KAUM); G. Shinohara, Y. Takata and M. Nakae (NSMT); Y. Sakurai (Okinawa Environmental Research Co. Ltd.); A. Murase and H.
Arihara (Tokyo University of Marine Science and Technology); M.
Mizukami (Fukuyama University); N. Akama, Y. Matsue, E. Nagao
and Y. Minagawa (professional fishermen on Hachijo-jima, Kozu,
Fukue, and Ishigaki-jima Island, respectively); Ogasawara Diving
Service KAIZIN; Ogasawara Flippers Inn; and the Fisheries Cooperative Associations on Chichi-jima Island, Minami-Izu-machi, Goto
Islands, Tanegashima Island and Yaku-shima Island. We also thankK. Kawabe (Tokyo Metropolitan Ogasawara Fisheries ResearchCenter) for useful information about ecology of E. fasciatus and T.
Yonezawa (Fudan University, China) for helping phylogenetic anddemographic analyses. We are grateful to G.S. Hardy (Ngunguru,
New Zealand) and E. Murdy (George Washington University) for
their help with English. This study was supported by a Grant-in-Aid
for Scientific Researches (A) and (B) from the Japan Society for the
Promotion of Science (A: 19208019, B: 24370041), a Grant-in-Aid
for Young Scientists (B) from the Ministry of Education, Science,Sports and Culture, Japan (22770089) and the Kuroshio Project of the
National Museum of Nature and Science.

References
Akaike H (1973) Information theory and an extension of the
maximum likelihood principle. In: Petrov BN, Csaki F (eds)Second international symposium on information theory. Akademiai Kiado, Budapest, Hungary, pp 267–281
Almany GR, Berumen ML, Thorrold SR, Planes S, Jones GP (2007)
Local replenishment of coral reef fish populations in a marine
reserve. Science 316:742–744
Avise JC (2000) The history and purview of phylogeography. In:
Avise JC (ed) Phylogeography: The history and formation of
species. Harvard University Press, Cambridge, pp 3–36
Bernardi G, Findley L, Rocha-Olivares A (2003) Vicariance anddispersal across Baja California in disjunct marine fish populations. Evolution 57:1599–1609
Bowen BW, Bass AL, Mus A, Carlin J, Robertson DR (2006a)
Phylogeography of two Atlantic squirrelfishes (Family Holocentridae): exploring links between pelagic larval duration andpopulation connectivity. Mar Biol 149:899–913
Bowen BW, Muss A, Rocha LA, Grant WS (2006b) Shallow mtDNA
coalescence in Atlantic pygmy angelfishes (Genus Centropyge)
indicates a recent invasion from the Indian Ocean. J Hered
97:1–12Briggs JC (1974) Chapter 8: Northern Hemisphere Warm-Temperate
Regions, and Chapter 9: Northern Hemisphere Cold-Temperate
Regions. In: Briggs JC (ed) Marine Zoogeography. McGrawHill, Inc., New York, pp 196–245, 247–297
Castelloe J, Templeton A (1994) Root probabilities for intraspecificgene trees under neutral coalescent theory. Mol Phylogenet Evol
3:102–113
Clement M, Posada D, Crandall KA (2000) TCS: A computer
program to estimate gene genealogies. Mol Ecol 9:1657–1659
Craig MT, Hastings PA (2007) A molecular phylogeny of the groupersof the subfamily Epinephelinae (Serranidae) with a revisedclassification of the Epinephelini. Ichthyol Res 54:1–17
Drummond AJ, Rambaut A (2007) BEAST: Bayesian evolutionary
analysis by sampling trees. BMC Evol Biol 7:214

Drummond AJ, Rambaut A, Shapiro B, Pybus OG (2005) Bayesiancoalescent inference of past population dynamics from molecular sequences. Mol Biol Evol 22:1185–1192
Dupanloup I, Schneider S, Excoffier L (2002) A simulated annealing
approach to define the genetic structure of populations. Mol Ecol
11:2571–2581
Eble JA, Toonen RJ, Bowen BW (2009) Endemism and dispersal:
comparative phylogeography of three surgeonfishes across the
Hawaiian Archipelago. Mar Biol 156:689–698
Excoffier L (2004) Patterns of DNA sequence diversity and genetic
structure after a range expansion: lessons from the infinite-islandmodel. Mol Ecol, 13:853–864
Excoffier L, Laval G, Schneider S (2005) ARLEQUIN ver. 3.0: an
integrated software package for population genetics data
analysis. Evol Bioinform Online 1:47–50
Excoffier L, Smouse PE, Quatro JM (1992) Analysis of molecularvariance inferred from metric distances among DNA haplotypes:
application to human mitochondrial DNA restriction data.
Genetics 131:479–491
Fu YX (1997) Statistical neutrality of mutations against population
growth, hitchhiking and background selection. Genetics
147:915–925
Fujikura K, Lindsay D, Kitazato H, Nishida S, Shirayama Y (2010)
Marine biodiversity in Japanese waters. PLoS One 5(8):e11836
Goloboff PA, Farris JS, Nixon KC (2008) TNT, a free program forphylogenetic analysis. Cladistics 24:774–786
Grant WS, Bowen BW (1998) Shallow population histories in deep
evolutionary lineages of marine fishes: insights from sardines
and anchovies and lessons for conservation. J Hered
89:415–426
Hasegawa M, Kishino H, Yano T (1985) Dating of the human-ape
splitting by a molecular clock of mitochondrial DNA. J Mol Evol
22:160–174
Herpending HC (1994) Signature of ancient population growth in a
low resolution mitochondrial DNA mismatch distribution. Hum
Biol 66:591–600
Hickey AJ, Lavery SD, Hannan DA, Scott-Baker C, Clements KD(2009) New Zealand triplefin fishes (family Tripterygiidae):contrasting population structure and mtDNA diversity within amarine species flock. Mol Ecol 18:680–696
Jobb G, von Haeseler A, Strimmer K (2004) Treefinder: a powerfulgraphical analysis environment for molecular phylogenetics.
BMC Evol Biol 4:18
Jobb G (2008) TREEFINDER Version October 2008. http://www.
treefinder.de/. Accessed August 2010
Jones GP, Milicich MJ, Emslie MJ, Lunow C (1999) Self-recruitment
in a coral reef fish population. Nature 402:802–803
Jordan DS (1901) The fish fauna of Japan, with observation on thegeographical distribution of fishes. Science, new series vol XIV
354:545–567
Kawabe K (2005) Embryonic development and effect of watertemperature on hatching of the Blacktip Grouper, Epinephelus
fasciatus. Aquaculture Sci 53:333–342
Kawabe K, Kato K, Kimura J (2000) Year round spawning of reared
blacktip grouper, Epinephelus fasciatus, in Chichi-jima, Ogasawara Islands, southern Japan. Aquaculture Sci 48:467–473
Kawabe K, Kato K, Kimura J, Saito M, Ando K, Kakiuchi K (1997)
Growth of reared blacktip grouper Epinephelus fasciatus in
Chichijima, Ogasawara Islands, Southern Japan. Aquaculture Sci
45:207–212
Kawabe K, Kohno H (2009) Morphological development of larva andjuvenile blacktip grouper, Epinephelus fasciatus. Fish Sci
75:1239–1251
Kimura M (2000) Paleogeography of the Ryukyu Islands. Tropics
10:5–24

123

K. Kuriiwa et al.
Kimura M (1980) A simple method for estimating evolutionary rates
of base substitutions through comparative studies of nucleotide
sequences. J Mol Evol 16:111–120
Kimura M, Weiss GH (1964) The stepping stone model of population
structure and the decrease of genetic correlation with distance.
Genetics 49:561–576
Kuriiwa K, Hanzawa N, Yoshino T, Kimura S, Nishida M (2007)
Phylogenetic relationships and natural hybridization in rabbitfishes (Teleostei: Siganidae) inferred from mitochondrial andnuclear DNA analyses. Mol Phylogenet Evol 45:69–80
Liu JX, Gao TX, Wu SF, Zhang YP, (2007) Pleistocene isolation in
the Northwestern Pacific marginal seas and limited dispersal in amarine fish, Chelon haematocheilus (Temminck and Schlegel,1845). Mol Ecol 16:275–288
Matsuura K, Senou H (2012) Introduction of Fishes in the Kuroshio
Current. In: Matsuura K (ed) Fishes in the Kuroshio Current.
Tokai University Press, Kanagawa, pp 3–16
Nakabo T (2002) Characteristics of the fish fauna of Japan andadjacent waters. In: Nakabo T (ed) Fishes of Japan with pictorial
keys to the species, English edition. Tokai University Press,
Tokyo, pp xliii–lii
Nishimura S (1992) General Remarks. In: Nishimura S (ed) Guide toseashore animals of Japan with color pictures and keys, vol. 1.
Hoikusha, Osaka, pp xi–xix
Ota H (1998) Geographic patterns of endemism and speciation in
amphibians and reptiles of the Ryukyu Archipelago, Japan, withspecial reference to their paleogeographical implications. Res
Popul Ecol 40:189–204
Palumbi SR (1994) Genetic divergence, reproductive isolation, andmarine speciation. Ann Rev Ecol Syst 25:547–572
Rambaut A, Drummond A (2009) TRACER v1.5. http://tree.bio.ed.
ac.uk/software/tracer/. Accessed January 2013
Randall JE, Heemstra PC (1991) Revision of Indo-Pacific groupers(Perciformes: Serranidae: Epinephelinae), with descriptions offive new species. Indo-Pacific Fishes No. 20. Bishop Museum,
Honolulu, Hawaii
Randall JE, Ida H, Kato K, Pyle RL, Earle JL (1997) Annotatedchecklist of the inshore fishes of the Ogasawara Islands. Natn Sci
Mus Monographs, No. 11. National Science Museum, Tokyo
Rice WR (1989) Analyzing tables of statistical test. Evolution
43:223–225
Riginos C, Nachman MW (2001) Population subdivision in marine
environments: the contributions of biogeography, geographicaldistance and discontinuous habitat to genetic differentiation in a
blennioid fish, Axoclinus nigricaudus. Mol Ecol 10:1439–1453
Rocha LA, Bass AL, Robertson R, Bowen BW (2002) Adult habitat
preferences, larval dispersal, and the comparative phylogeography of three Atlantic surgeonfishes (Teleostei: Acanthuridae).
Mol Ecol 11:243–252
Rocha LA, Rocha CR, Robertson DR, Bowen BW (2008) Comparative phylogeography of Atlantic reef fishes indicates both origin
and accumulation of diversity in the Caribbean. BMC Evol Biol
8:157
Rogers AR (1995) Genetic evidence for a Pliocene population
explosion. Evolution 49:608–615
Rogers AR, Harpending H (1992) Population growth makes waves inthe distribution of pairwise genetic differences. Mol Biol Evol
9:552–569
Ronquist F, Teslenko M, Van Der Mark P, Ayres DL, Darling A,
Höhna S, Larget B, Liu L, Suchard MA, Huelsenbeck JP (2012)
MrBayes 3.2: Efficient Bayesian phylogenetic inference andmodel choice across a large model space. Syst Biol 61:539–542
Schneider S, Excoffier L (1999) Estimation of past demographic
parameters from the distribution of pairwise differences when123

the mutation rates vary among sites: application to humanmitochondrial DNA. Genetics 152:1079–1089
Senou H (2002) Family Serranidae. In: Nakabo T (ed) Fishes of Japan
with pictorial keys to the species, English edition. Tokai
University Press, Tokyo, pp 690–731
Senou H, Shinohara G, Matsuura K, Furuse K, Kato S, Kikuchi T(2002) Fishes of Hachijo-jima Island, Izu Islands Group, Tokyo,
Japan. Mem Nat Sci Mus Tokyo 38:195–237
Senou H, Morita Y, Morishita O (2003) Notes on the distribution of asurgeonfish Ctenochaetus hawaiiensis (Perciformes: Acanthuridae). I.O.P. Diving News 14(10):2–4
Senou H, Matsuura K, Shinohara G (2006) Checklist of fishes in the
Sagami Sea with zoogeographical comments on shallow water
fishes occurring along the coastlines under the influence of theKuroshio Current. Mem Nat Sci Mus Tokyo 41:389–542
Shimada K (2002) Family Chaetodontidae. In: Nakabo T (ed) Fishesof Japan with pictorial keys to the species, English edition. Tokai
University Press, Tokyo, pp 884–897
Shirai S, Kuranaga R, Sugiyama H, Higuchi M (2006) Population
structure of the sailfin sandfish, Arctoscopus japonicas (Trichondontidae), in the Sea of Japan. Ichthyol Res 53:357–368
Springer VG (1982) Pacific Plate biogeography, with special reference to shorefishes. Smithson Contrib Zool 367. Smithsonian
Institution Press, Washington.
Swearer SE, Caselle JE, Lea DW, Warner RR (1999) Larval retention
and recruitment in an island population of a coral-reef fish.
Nature 402:799–802
Tajima F (1989) Statistical method for testing the neutral mutation
hypothesis by DNA polymorphism. Genetics 123:585–595
Tamura K, Nei M (1993) Estimation of the number of nucleotide
substitutions in the control region of mitochondrial DNA in
humans and chimpanzees. Mol Biol Evol 10:512–526
Tamura K, Peterson D, Stecher G, Nei M, Kumar S (2011) MEGA5:
Molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol 28:2731–2739
Tavaré S (1986) Some probabilistic and statistical problems in theanalysis of DNA sequences. Lect Mathemat Life Sci (AmMathemat Soc) 17:57–86
Tanabe AS (2007) KAKUSAN: a computer program to automate theselection of a nucleotide substitution model and the configurationof a mixed model on multilocus data. Mol Ecol Notes 7:962–964
Tanabe AS (2008) Phylogears Version 1.5.2010.03.24. http://www.
fifthdimension.jp/. Accessed August 2010
Tanaka S (1931) On the distribution of fishes in Japanese water. J Fac
Sci, Imp Univ Tokyo Sec IV Zool 3(1):1–90
Templeton AR, Crandall KA, Sing CF (1992) A cladistic-analysis ofphenotypic associations with haplotypes inferred from restriction
endonuclease mapping and DNA-sequence data. 3. Cladogram
estimation. Genetics 132:619–633
Ujiié H, Ujiié Y (1999) Late Quaternary course changes of theKuroshio Current in the Ryukyu Arc region, northwestern Pacific
Ocean. Mar Micropaleontol 37:23–40
Ujiié Y, Ujiié H, Taira A, Nakamura T, Oguri K (2003) Spatial andtemporal variability of surface water in the Kuroshio source
region, Pacific Ocean, over the past 21,000 years: evidence fromplanktonic forminifera. Mar Micropaleontol 49:335–364
Vos RA (2003) Accelerated likelihood surface exploration: thelikelihood ratchet. Syst Biol 52:368–373
Yang Z (2007) PAML 4: phylogenetic analysis by maximum
likelihood. Mol Biol Evol 24:1586–1591
Yoshigou H (2004) The fishes found from tide-pools and surge
sublittral zone in the Minami-daito (South-Borodino) Island,
Japan. Misc Rep Hiwa Mus Natl Hist Hiroshima 43:1–51