Genetic patterns in peripheral marine populations of the fusilier fish Caesio cuning within the Kuroshio Current

Abstract Aim Mayr's central‐peripheral population model (CCPM) describes the marked differences between central and peripheral populations in genetic diversity, gene flow, and census size. When isolation leads to genetic divergence, these peripheral populations have high evolutionary value and can influence biogeographic patterns. In tropical marine species with pelagic larvae, powerful western‐boundary currents have great potential to shape the genetic characteristics of peripheral populations at latitudinal extremes. We tested for the genetic patterns expected by the CCPM in peripheral populations that are located within the Kuroshio Current for the Indo‐Pacific reef fish, Caesio cuning. Methods We used a panel of 2,677 SNPs generated from restriction site‐associated DNA (RAD) sequencing to investigate genetic diversity, relatedness, effective population size, and spatial patterns of population connectivity from central to peripheral populations of C. cuning along the Kuroshio Current. Results Principal component and cluster analyses indicated a genetically distinct lineage at the periphery of the C. cuning species range and examination of SNPs putatively under divergent selection suggested potential for local adaptation in this region. We found signatures of isolation‐by‐distance and significant genetic differences between nearly all sites. Sites closest to the periphery exhibited increased within‐population relatedness and decreased effective population size. Main Conclusions Despite the potential for homogenizing gene flow along the Kuroshio Current, peripheral populations in C. cuning conform to the predictions of the CCPM. While oceanography, habitat availability, and dispersal ability are all likely to shape the patterns found in C. cuning across this central‐peripheral junction, the impacts of genetic drift and natural selection in increasing smaller peripheral populations appear to be probable influences on the lineage divergence found in the Ryukyu Islands.


| INTRODUC TI ON
Within the spatial distribution of a species, peripheral populations can be prone to edge effects that significantly alter their genetic characteristics relative to central counterparts. This phenomenon is summarized in the central-peripheral population model (CPPM; Mayr, 1963). Under the CPPM, populations at the center of a species range are contiguous, abundance is regulated by density-dependent factors, and gene flow is multidirectional.
These combined forces result in larger central populations with high genetic diversity. In contrast, peripheral populations are often fragmented and gene flow occurs from a single direction.
Empirical research indicates that the majority of species approximately conform to the CPPM but studies of marine taxa lag behind those of terrestrial taxa (Ecker, Samis, & Lougheed, 2008). Most marine organisms disperse via a larval pelagic phase, and a large range of dispersal potential exists among marine species due to differences in duration of the pelagic period, ability to orient, and oceanographic patterns complicating the predictability of population connectivity (Cowen, Gawarkiewicz, Pineda, Thorrold, & Werner, 2007;Kinlan, Gaines, & Lester, 2005;Selkoe, Gaggiotti, Bowen, & Toonen, 2014;Shanks, 2009;Weersing & Toonen, 2009). Scientists consider species with disjunct peripheral populations and low genetic diversity to have higher evolutionary value than those with continuous peripheral populations (Bunnell, Campbell, & Squires, 2004;Lesica & Allendorf, 1995). Therefore, it is important to increase our understanding of marine species and geographic regions prone to edge effects.
One region of interest for marine species is the Ryukyu Islands, which mark the northern periphery of many ubiquitous coral reef species' ranges in the western Pacific Ocean (Randall, 1998;Veron & Minchin, 1992). This trend is generally attributed to the steady decrease in sea surface temperature to beyond the thermodynamic threshold of tropical reef organisms at increasingly higher latitudes (Crossland, 1984;Dana, 1843;Jokiel & Coles, 1977;Rosen, 1975;Vaughan, 1919). Hermatypic corals are obligatory habitat for many tropical reef species and their prey, and at the latitude of the southernmost main island of Japan, shallow water reefs are classified as marginal due to low surface temperature, low aragonite saturation states, and low light (Kleypas, McManus, & Meñez, 1999). Without supportive habitat and beyond the extent of their fundamental niche, most tropical reef organisms do not succeed in establishing populations above the Ryukyus, though range expansions are beginning to be documented with rising global sea surface temperatures (Yamano, Sugihara, & Nomura, 2011).
The steady, northerly stream of warm surface water flowing through the Ryukyu Islands that results in the persistence of coral reef organisms at higher-than-normal latitudes is provided by a powerful western-boundary current known as the Kuroshio Current (Veron, 1992;Yamano, Hori, Yamauchi, Yamagawa, & Ohmura, 2001), but the role that western boundary currents such as the Kuroshio Current and its Atlantic Ocean counterpart, the Gulf Stream, may play in shaping the genetic patterns in the northern peripheral populations of tropical species relative to central counterparts remains mostly unexamined. The only taxon sampled across the flow path of the Kuroshio Current to address questions about genetic connectivity and diversity in central and peripheral populations is passive propagule-dispersing seagrass (Arriesgado et al., 2015(Arriesgado et al., , 2016Kurokochi et al., 2016;Nakajima et al., 2014) leaving a majority of the tropical species and ecosystems whose peripheral populations are found in the Ryukyu Islands unstudied.
The pathway of the Kuroshio Current is an intriguing setting for examining central-peripheral patterns in shallow water tropical marine populations since it contains both the topography to support disjunct populations and the oceanography to support continuous populations. The Kuroshio is formed when the Northern Equatorial Current (NEC) hits and bifurcates north and south at ~13°N along the eastern coastline of the Philippines ( Figure 1; Nitani, 1972;Toole, Millard, Wang, & Pu, 1990). It flows past nearly 1,000 km of contiguous shoreline habitat and crosses a large expanse of deep water to the disjunct, sub-tropical island reef systems that make up the Ryukyu Islands before finally turning east off the coast of Honshu, Japan to form the North Pacific Current. However, it can take as little as one effective migrant per generation to nullify disruptive genetic drift between two populations (Spieth, 1974), and since the Kuroshio Current can reach mean maximum surface velocities of ~1.2 m/s (~104 km/day) (Yang et al., 2015) there is great potential for genetic continuity via the dispersive larvae characteristic of most marine organisms.
As one of the many tropical species whose northern distribution ends in the Ryukyu Islands, the Redbelly Yellowtail Fusilier, Caesio cuning (Bloch, 1791), has the potential for either disjunct or continuous peripheral populations within the Kuroshio Current.
Like the majority of reef organisms, C. cuning has a bipartite life history beginning as pelagic larvae and settling on coral reefs as juveniles. In contrast to other fusiliers, such as those in the genera Pterocaesio and Dipterygonotus, it has been noted that C. cuning larvae are primarily "coastal", being nearly always found over the mid to inner continental shelf (Reader & Leis, 1996).
Observations of eleven C. cuning larvae in situ also showed some evidence that larvae can detect and orient themselves toward reefs (Leis & Carson-Ewart, 2003). Adults are non-migratory and dependent on reef structure for protection at night. These observations suggest that long distance dispersal in C. cuning is unlikely without the presence of a strong oceanographic conduit such as the Kuroshio Current. And as an exploited fishery throughout its range, a better understanding of whether C. cuning peripheral populations are continuous or disjunct will be useful for effective management.

| MATERIAL S AND ME THODS
Pectoral fin and muscle tissue were collected from 307 fish at regional markets and landings from three sites below the latitude which marks the standard tropical-subtropical boundary (23.5°N) along the east coast of the Philippines and two sites at the northern edge of the species range above 23.5°N in the Ryukyu Islands of Japan ( Figure 1, Table 1). Tissues were stored in 95% ETOH or DNA/RNA Shield™ (Zymo Research), and DNA was extracted using E-Z 96 ® Tissue DNA Kits (Omega Bio-tek, Inc.). Elutions containing high-weight DNA were quantified with a fluorescence microplate reader, and aliquots containing 100 ng of DNA were cleaned using AMPureXP beads (Beckman Coulter, Inc.) for input into library preparation.

| RAD library preparation and sequencing
Restriction site-associated DNA (RAD) libraries were prepared following a modified ezRAD protocol (Toonen et al., 2013). Genomic DNA was cut at 5'-GATC sites using the isoschizomers MboI and Sau3AI (New England Biolabs). Cleaned template DNA bound by AMPureXP beads was eluted in 21.5 µl of water and digested overnight with 2.5 uL of CutSmart ® Buffer and 2.5 U of each restriction enzyme in 25 µl reactions. Digested DNA was rebound to beads using a 2× 3 M NaCl, 20% PEG solution (Faircloth & Glenn, 2014;Fisher et al., 2011) and cleaned before being input into an Illumina

| Read processing, SNP filtering, and outlier detection
Forward and reverse reads were truncated to 90 bp and quality filtered using the Stacks subprogram process_radtags (Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013) before being input into the dDocent pipeline (Puritz, Hollenbeck, & Gold, 2014) for de novo assembly, read mapping, and variant calling. Default dDocent parameters were modified based on performance, including increasing the CD-HIT sequence similarity parameter (-c) to 0.92 and increasing the minimum base phred score (-q) to 20 for an allele to be included in variant calling with Freebayes v1.0.2-58-g054b257 in order to decrease the overall number of RAD contigs that had more than two haplotypes present after mapping.
Filtering parameters were modified from recommendations in dDocent content (https://github.com/jpuritz/dDocent/). Preliminary filtering to remove samples with high amounts of missing data and erroneous calls was done with VCFtools (Danecek et al., 2011).
Conditions under which variants were removed were as follows: quality value lower than 30, genotyped in fewer than 95% of individuals, a minor allele frequency less than 0.05, a mean coverage depth less than five, and sites with more than two alleles. Additional filtering with the program vcflib (https://github.com/ekg/vcflib) removed loci with a heterozygote allele balance <0.25 and>0.75, reads from both strands, large variation in mapping quality among alleles, an alternate allele only supported by unpaired reads, and a mean depth of well above the majority distribution of depths across loci since high coverage can lead to inflated quality scores or false heterozygotes (Li, 2014). Remaining variants were decomposed to remove indels. BayeScan uses a Bayesian approach to selection detection that implements the multinomial-Dirichlet model and was run with default parameters. A false discovery rate (FDR) correction (α = 0.05) was set in both programs, and any loci identified as outliers by either program after correction were removed to produce a neutral panel of SNPs for our analyses. Both forward and reverse tags associated with outlier loci were run through the megablast search tool in BLASTN v2.6.1 to examine similarity to NCBI archived sequences in the nr/nt and wgs databases (Morgulis et al., 2008;Zhang, Schwartz, Wagner, & Miller, 2000).

| Genetic differentiation and spatial clustering
Multiple approaches were used to examine the level of genetic connectivity among sites. Pairwise genetic differentiation between populations (F ST ) was generated from putatively neutral SNPs in the program GenoDive v2.0b27 (Meirmans & Van Tienderen, 2004).
Significance of F ST values was tested using 10,000 permutations, and p-values were corrected using Benjamini & Hochberg's method of FDR correction (Benjamini & Hochberg, 1995). Genetic variability within and among sites was examined in both neutral and putatively adaptive loci with a principal components analysis (PCA) in the R package "adegenet" (Jombart, 2008;Jombart & Ahmed, 2011). All remaining analyses were run exclusively on the panel of neutral loci.
The power of different methods to estimate the number of discrete genetic populations within a dataset can vary based on a species' dispersal characteristics and the spatial distribution of sampling (Murphy, Evans, Cushman, & Storfer, 2008;Schwartz & McKelvey, 2008), therefore the use of multiple methods to detect barriers is recommended (Blair et al., 2012). The program STRUCTURE v2.3.4 (Falush, Stephens, & Pritchard, 2003;Pritchard, Stephens, & Donnelly, 2000) was used to identify distinct genetic clusters and assign individuals to populations under the assumption of population admixture and correlated allele frequency. Runs consisted of a burn-in period of 50,000 MCMC iterations followed by 100,000 iterations for inferred K of 1 to 5 and replicated 10 times for consistency. Population identity was not used as an a priori parameter for clustering. Results were collated in STRUCTURE HARVESTER (Earl & vonHoldt, 2012), and the log probability and ΔK statistic (Evanno, Regnaut, & Goudet, 2005) were used to determine the most likely number of clusters. Replicate Q-matrices for the optimal K were processed with the Greedy algorithm in CLUMPP v1.1.2 (Jakobsson & Rosenberg, 2007)  Samples from each location shared the same spatial coordinates so to allow individuals with the same coordinates to be assigned to different populations, UTM coordinate uncertainty was set to 1 km.
The modal K from these initial runs was fixed for 20 additional replicates, and the run with the highest posterior probability was used to assign individuals to clusters.
Distinguishing spatial clustering from clinal genetic differentiation can be challenging since tests for one can be biased by the presence of the other (Guillot, Leblois, Coulon, & Frantz, 2009;Meirmans, 2012;Schwartz & McKelvey, 2008). One suggested method of distinguishing a pattern of isolation-by-distance (IBD) from spatial clustering is running partial Mantel tests. Shortest overwater distances between collection sites were calculated in ArcGIS v10.1 using the Cylindrical Equal Area projected coordinate system and a 500 m cell size (2,500 m 2 ). Distances were estimated using the cost distance tool, assigning an equal "cost" to each cell of water

| Genetic diversity, effective population size, and relatedness
To examine sites for the presence of edge effects, we assessed genetic diversity, effective population size, and relatedness within Coefficients of relatedness (r) for individuals at each site were estimated using a moment estimator and two maximum-likelihood estimators in the R package "related" (Pew, Muir, Wang, & Frasier, 2015), which have been shown to be robust methods for examining relatedness (Attard, Beheregaray, & Möller, 2018;Wang, 2011Wang, , 2014. One pair of individuals, San_007 and San_008, had a r value of 0.97 which is consistent with the same individual being sampled twice so San_008 was removed from all further analyses (Supporting Information, Table S1).   the fish from Ishigaki were also distinct, resulting in three clusters-one comprised of the samples from Okinawa, one comprised of the samples from Ishigaki, and one comprised of the Philippine sites ( Figure 2b). PC one separates the two Ryukyu Island sites from the Philippines sites and accounts for 11.9% of the variance in these SNPs. No signal of lane or sequencer effects was found in any PCs.

| Genetic differentiation and spatial clustering
STRUCTURE results were consistent with those from the PCAs, with all but one fish from the peripheral population of Okinawa belonging to an exclusive genetic cluster (Figure 3). Both the log-likelihood [L(K)] and ΔK indicated that K = 2. The remaining Okinawan fish exhibited a profile similar to those from Ishigaki, which were assigned to a second cluster with all individuals from the central Philippines sites (Figure 3). The GENELAND run with the largest log-likelihood score assigned the individuals from the three Philippine sites to one cluster, the individuals from Ishigaki to a second cluster, and the individuals from Okinawa to a third cluster.
A basic Mantel test was significant for correlation between genetic and geographic distance (r = 0.6061, p = 0.0083, Figure 4)

| Genetic diversity, effective population size and relatedness
Summary statistics by site are reported in Table 1 Estimates of effective population size (N e ) within sites ranged from 2,258-5,731 in Okinawa (peripheral) to "infinite" in Atimonan and Guiuan (central , Table 4). An "infinite" estimate is often an indication that the mean sample size (S) is too small to generate an estimate of N e using the LD N e method (Do et al., 2014;Waples & Do, 2010). When all pairwise locus comparisons were included in estimates (r 2 < 1), the smallest N e was measured in Okinawa, with increasing estimates of N e in Ishigaki and Santa Ana. Given that sample sizes were similar among all sites (S = 30.4-38.5), and these were adequate to generate finite estimates of N e in Okinawa, Ishigaki, and Santa Ana but not Atimonan and Guiuan, it indicates that the true N e in these two southern sites is larger than the true N e of the northern sites, which conforms to our expectations. Exclusion of pairwise locus comparisons using a cutoff of r 2 ≤ 0.5 did affect estimates of N e , indicating the presence of some linkage between loci.
The peripheral population of Okinawa had the most genetically re-  Ritland, 1996) and one pair of individuals from Okinawa r = 0.159, Milligan, 2003, r = 0.158, Wang, 2007, and r = 0.156, Ritland, 1996 was approximately that expected from first cousins. An additional 11 Okinawa-Okinawa pairs, an Okinawa-Ishigaki pair, and an Atimonan-Guiuan pair had relatedness estimates with an upper 95% confidence limit at or above that expected from first cousins. The related pairs found within Ishigaki and Okinawa contributed to an increase in average relatedness measured among individuals within sites in the Ryukyu Islands, with the largest average relatedness measured within Okinawa ( Figure 5). Mean coefficients of relatedness for all estimators are reported in Table S3 in Supporting Information.   authors note that ghost populations in simulated data are extremely rare and suggest the presence of them can be a clue that data depart from model assumptions . A departure such as isolation-by-distance, which partial mantel tests indicated is present across the five sampled sites, could result in overestimating the number of genetic clusters.

| Genetic differentiation within the Kuroshio Current
The presence of isolation-by-distance between sites supports evidence from observations of larval behavior that C. cuning do not regularly disperse long distances. A previous study using mitochondrial DNA found no evidence for IBD though long-distance dispersal was hypothesized to be rare from the low levels of mixing observed between divergent mitochondrial clades (Ackiss et al., 2013). This is the first time IBD has been documented in this species and is likely due to the resolution provided by a large SNP dataset. Notes. Effective population size estimates (N e ) and confidence intervals (CI) were generated in R using the r 2 values from the linkage disequilibrium method in NeEstimator with a minor allele frequency cutoff of 0.05. N e for each location or combination of locations was calculated with all pairwise comparisons (r 2 < 1) and using only comparisons with r 2 ≤ 0.5. The NEC bifurcation occurs between the latitudes of 11°-16.5°N (Qiu & Lukas, 1996)  The results of this study indicate that genetic patterns from central to peripheral populations of C. cuning within the Kuroshio Current are influenced by both dispersal ability and local adaptation in marginal habitat. With the predominance of this western-boundary current and the shelf topography in this region, no known physical barrier has limited pelagic larval dispersal, even during the low sea levels of the Pleistocene epoch (Voris, 2000).

| Edge effects in peripheral populations
Given the dominance of a single lineage in the Philippines, continuous habitat appears to be influential enough to prevent large amounts of genetic divergence. Once this continuous shoreline F I G U R E 5 Relatedness among C. cuning within sites. Mean coefficients of relatedness (r) generated with three different estimators in the R package "related." Standard error of the mean is denoted by bars  Ritland (1996) Milligan (2003) Wang (2007) Site habitat gives way to disjunct islands, patterns of isolation-by-distance, and peripheral divergence arise. Enough dispersal from the Philippines to Ishigaki appears to be occurring to make individuals in this island more similar to those in the Philippines than Okinawa, even though Okinawa is nearly 200 km closer to Ishigaki. However, what roles successively decreasing effective population sizes across these disjunct stepping-stones and localized selective pressures in increasingly subtropical habitat play in shaping this pattern cannot be determined from these data.
Studies of additional taxa with a range of dispersal potential and peripheral genetic diversity or comparative studies examining differences in gene expression in individuals in central and peripheral populations will be needed to further untangle the relative influences of dispersal ability and selective pressure in shaping patterns in peripheral populations of tropical species in the Kuroshio Current.

ACK N OWLED G M ENTS
We are grateful to A. Lizano for his persistence collecting from

CO N FLI C T O F I NTE R E S T
None declared.