Observations of Migrant Exchange and Mixing in a Coral Reef Fish Metapopulation Link Scales of Marine Population Connectivity

Much progress has been made toward understanding marine metapopulation dynamics, largely because of multilocus microsatellite surveys able to connect related individuals within the metapopulation. However, most studies are focused on small spatial scales, tens of kilometers, while demographic exchange at larger spatial scales remains poorly documented. Additionally, many small-scale demographic studies conflict with broad-scale phylogeographic patterns concerning levels of marine population connectivity, highlighting a need for data on more intermediate scales. Here, we investigated demographic recruitment processes of a commercially important coral reef fish, the bluespine unicornfish ( Naso unicornis ) using a suite of mitochondrial DNA (mtDNA) and microsatellite markers. Sampling for this study ranged across the southern Marianas Islands, a linear distance of 250 km and included 386 newly settled postlarval recruits. In contrast with other studies, we report that cohorts of recruits were genetically homogeneous in space and time, with no evidence of temporally stochastic sweepstakes reproduction. The genetic diversity of recruits was high and commensurate with that of the adult population. In addition, there is substantial evidence that 2 recruits, separated by 250 km, were full siblings. This is the largest direct observa-tion of dispersal to date for a coral reef fish. All indications suggest that subpopulations of N. unicornis experience high levels of demographic migrant exchange and metapopulation mixing on a spatial scale of hundreds of kilometers, consistent with high levels of broad-scale genetic connectivity previously reported in this species.

Benthic marine organisms with a pelagic reproductive strategy, such as many coral reef fishes, are subject to the highly advective environment of the ocean. After spawning, the buoyant eggs and pelagic-dwelling larvae are transported away from their natal, shallow waters by oceanographic forces (Domeier and Colin 1997;Pringle and Wares 2007;Siegel et al. 2008;Cowen and Sponaugle 2009). Formerly, it was believed that these oceanographic processes led to passive dispersal and spatial disassociation of larval propagules, resulting in randomly mixed larval pools. The reproductive outputs of local adult populations, therefore, appeared decoupled from rates of larval recruitment and populations were deemed "demographically open," such that populations were predicted to randomly exchange migrants on a massive scale through ocean circulations (Doherty 1991;Sale 1991).
Research has since shown, however, that most marine populations do not conform to traditional notions of demographic openness and that events taking place during the pelagic larval phase tend to result in the genetic heterogeneity of recruitment cohorts, rather than a homogenized larval pool. For instance, many studies report that a significant portion of pelagic larvae self-recruit to their natal habitats (Jones Horne et al. • Migrant Exchange in a Coral Reef Fish Metapopulation et al. 1999Swearer et al. 1999;Almany et al. 2007;Gerlach et al. 2007;Patterson and Swearer 2007;Christie et al. 2010a;Small and Wares 2010;Horne et al. 2011;Buston et al. 2012;Priest et al. 2012), refuting the notion that populations are seeded primarily from external sources. Additionally, larval cohorts sampled from the same site exhibit temporal genetic structure, referred to as chaotic genetic patchiness (Hellberg et al. 2002), attributed to stochastic oceanographic forces (Selkoe et al. 2006;Pringle and Wares 2007;Siegel et al. 2008;Cowen and Sponaugle 2009), kinship aggregations (Avise and Shapiro 1986;Selkoe et al. 2006;Buston et al. 2009;Bernardi et al. 2012), and temporal stochastic variation in reproductive success (Sweepstakes reproduction- Hedgecock et al. 2007). A temporal pattern, such as sweepstakes reproduction, where each cohort of recruits is stochastically chosen from a small subset of potential parents located somewhere in the metapopulation, argues that migrant exchange between any 2 subpopulations is likely to be irregular and potentially limiting in terms of demographic impact. Therefore, marine populations that experience self-recruitment and sweepstakes reproduction may, to a large degree, be demographically closed. However, the extent to which marine populations exchange demographically relevant amounts of migrants at various spatial scales continues to be a subject of debate.
Through direct genetic observations of parent-offspring pairs, it is now known that dispersal in coral reef fishes is demographically relevant at scales of tens to more than a hundred kilometers Christie et al. 2010b;Saenz-Agudelo et al. 2011;Harrison et al. 2012). The rapid spread of invasive exotic reef fishes, such as lionfish (Pterois volitans) in the Caribbean (Betancur-R. et al. 2011) and the snapper Lutjanus kasmira in the Hawaiian archipelago (Gaither et al. 2010), supports the idea that substantial demographic exchange can occur across larger distances also (Mora et al. 2012). However, with only a few exceptions (see Hepburn et al. 2009;Christie et al. 2010b;Jones et al. 2010), studies to date have not attempted to detect relatedness in individuals at distances greater than a few kilometers, nor have they examined how genetic composition varies over time. Knowledge of the composition of pools of dispersing individuals is fundamental to the understanding of metapopulation dynamics (Slatkin 1977), which in turn is fundamentally important for spatially explicit conservation initiatives, such as networks of marine protected areas (Hastings and Botsford 2006;Kritzer and Sale 2006;Jones et al. 2009). In this regard, empirical data on demographically relevant genetic patterns, at all spatial scales, is desired (Botsford et al. 2001;Fogarty and Botsford 2007).
A useful coral reef fish species for studying recruitment processes is the acanthurid Naso unicornis, the bluespine unicornfish. The settlement-stage larvae of this species are large (50-70 mm), easy to identify and settle specifically in shallow reef-flat lagoons, occasionally in large discrete pulses that allow for robust sampling. For these attributes, N. unicornis has previously been the target of demographic studies of coral reef fishes Doherty et al. 2004) and is 1 of a small number of widespread species that have been genetically surveyed across the Indo-Pacific region (Horne et al. 2008). What is already known about N. unicornis indicates that it has a long pelagic larval duration (approximately 71 days- Wilson and McCormick 1999), with a right-skewed distribution in the age of pelagic larvae, suggesting that it can delay settlement until after 95 days (Doherty et al. 2004), and that it experiences high levels of long-term gene flow across its species distribution (Horne et al. 2008; see also Horne and van Herwerden 2013). Therefore, this species is a good candidate for investigating demographic exchange in coral reef fishes beyond what other studies have heretofore attempted.
A study by  used 19 allozyme loci and compared the genetic composition of recruiting N. unicornis larvae to adults and postsettlement juveniles on Moorea and discovered temporal stochasticity in this species. This outcome suggests that even in N. unicornis, which is presumed to be among the most dispersive of coral reef fishes, demographically significant migrant exchange at anything greater than local spatial scales may be restricted. The large recruitment pulses observed in this species might, therefore, be the offspring of a few highly fecund progenitors  favored by sweepstakes events rather than a mixed larval pool. Yet, at the same time, these small-scale demographic patterns do not necessarily preclude high levels of population connectivity (Armsworth 2002;Hepburn et al. 2009). Thus, whether there is a mismatch between large-scale evolutionary patterns of gene flow and ecological patterns of demographic exchange at smaller scales, in this and other perceptually dispersive marine organisms, remains an unresolved but important research priority (Hixon 2011).
In the present study, newly recruited postlarval juveniles (recruits) of N. unicornis were collected from 6 inner reef flat sites across 2 islands in the southern Marianas archipelago and included 4 major recruitment pulses in 2008 and 2009. Genetic variation was surveyed for 1 mtDNA marker and 12 microsatellite loci, which are more polymorphic and more sensitive to ecological patterns than allozymes (Waples and Gaggiotti 2006). Furthermore, in addition to other motivations for studying N. unicornis, this is a fishery species in many parts of the Indo-Pacific region and makes up a sizeable portion of total fishery yields (Rhodes et al. 2008;Houk et al. 2012). It also plays a key functional role in coral reef ecosystems-as a consumer of coral-suppressing macroalgae that are unpalatable to most marine herbivores (Hoey and Bellwood 2009). Therefore, demographic investigations of N. unicornis provide data relevant to fisheries and conservation management of this species.

Sampling
We sampled 386 N. unicornis recruits from reef flat habitats at 5 sites around Guam: 2 from the west coast (Tanguison and Asan), 3 from the east coast (Cocos Lagoon, Ipan Beach, and Pago Bay), and 1 site on western Saipan (Figure 1). Our sampling coincided with 4 major recruitment pulses that occurred in June 2008, July 2008, September 2008, and September 2009. The first pulse, in particular, was quite large. During a similar event in French Polynesia, involving the same species, Doherty et al. (2004) reported that as many as 10 000 recruits per square kilometer can settle on the reef per night at the height of settlement. Nevertheless, postsettlement mortality removes about 80% of these individuals in the first few days (Doherty et al. 2004). All sites on Guam were sampled at least once per year in 2008 and 2009 (Tables 1 and 4). Samples were gathered by hand, at night, by snorkel at depths no greater than 1 m. Immediately following collection, samples were transported to the University of Guam Marine Laboratory, where a small fin clip was taken and preserved in 80% ethanol. Samples were weighed and measured for fork length. During the same period, 176 adult fish were collected from around Guam either by spearing or purchase from local commercial vendors. Adults were treated as a single population, considering that spawning events are likely to draw individuals from large areas (see Discussion).

Laboratory Procedures
Total genomic DNA was extracted using a Chelex extraction method (Walsh et al. 1991). First, a segment of the mitochondrial control region (mtCR) was amplified using the genusspecific NA1 primers of Klanten et al. (2007). Polymerase chain reaction (PCR) parameters and thermocycling profiles for the mtCR are identical to those described in the study by Horne et al. (2008). PCR products were purified and sequenced with the NA1 forward primer (5' AGC ATT CTG AAC TAA ACT AC 3') at Macrogen Sequencing Service, Seoul, South Korea. MtCR haplotype alignments from this study are publicly available from Dryad (doi:10.5061/dryad. hm796). Sequences that were difficult to read were not included in this study.
Second, all samples were genotyped at 12 microsatellite loci, Nuni 01-Nuni 12, described by Horne et al. (2010). PCR amplification conditions of the microsatellite loci were identical to those described by Horne et al. (2010). PCR products were purified using ethanol and ammonium acetate precipitation and read using Amersham MegaBACE instrumentation at the James Cook University Genetic Analysis Facility. Seventy-seven samples (14%), including redos, were genotyped twice for quality control, with consistent results. Microsatellite scores from this study are available from Dryad (doi:10.5061/dryad.hm796).

Genetic Analyses of Microsatellite Loci
Microsatellite peaks were scored using the program FRAGMENT PROFILER version 1.2 (Amersham Biosciences). Exact tests for departure from Hardy-Weinberg equilibrium (HWE) were conducted in GENEPOP version 4.0.10 (Rousset 2008) and linkage disequilibrium (LD) was mapped using the Markov chain algorithm, a dememorization of 10 000, with 20 batches and 5000 iterations per batch. Loci that did not conform to HWE may be under selection and were not used for downstream analyses, but this was based solely on allele frequencies from the single adult population on Guam because recruitment cohorts are not reproductive populations, potentially contain kinship aggregations, and may display allele frequencies that are otherwise misleading in this regard. The presence of null alleles, large allele dropout, stuttering, and other genotyping errors were assessed in MICROCHECKER version 2.2.3 (van Oosterhout et al. 2004). Number of alleles, allelic richness, number of private alleles, and the observed and expected heterozygosities based on Hardy-Weinberg proportions were estimated in FSAT version 2.9.3 (Goudet 2001) and in GENALEX version 6.5 (Peakall and Smouse 2012). Pairwise standardized D values (Jost 2008) were calculated for each site in GENALEX with 10 000 permutations of the data. To safeguard against type I error in pairwise comparisons, we applied the false discovery rate of Benjamini and Yekutieli (2001), as recommended by Narum (2006).

Discriminant Analysis of Principal Components
Multivariate discriminant analysis of principal components (DAPC) scatter plots were created in R (R Development Core Team 2009; http://www.Rproject.org), using the package adegenet (Jombart 2008;Jombart et al. 2010), to visually represent genetic patterns among N. unicornis recruits and to give the best genetic discrimination of these groups as they were collected in the field. However, due to a lack of spatial or temporal genetic structure (see Results), there was no justification for assigning individuals to predefined sample locations. Instead the "find.clusters" function from the R package adegenet (Jombart et al. 2010) was used to detect genetic clustering without prior group information. This analysis is equivalent to that done more conventionally in the program STRUCTURE (Pritchard et al. 2000) but requires much less computational time. First, data were transformed into unrelated variables using principal components analysis (PCA). For this analysis, data were scaled and all missing data were assigned to the mean of PCA. Next, a number of PCs were retained as predictors for discriminant analysis. There are no strict guidelines for determining how many PCs should be retained during this dimensions-reduction step, but it is a compromise between the statistical power of more PCs and the stability of results (Jombart et al. 2010 and references therein). For the purposes of this study, 100 PCs were retained, containing 80% of the variation of the data. The correct number of demes (k) was selected based on likelihood score and the Bayesian information criterion with 10 000 iterations.

Relatedness of Recruitment Cohorts
To investigate sibling relationships, pairwise relatedness coefficients were calculated for all recruits. Three different estimators were used: Queller and Goodnight's (1989) coefficient of relatedness, Lynch and Ritland's (1999) estimator, and the estimator of Wang (2002). There is no clear consensus about the estimator that most accurately identifies siblings under any given demographic scenario (Oliehoek et al. 2006;Konovalov and Heg 2008), but all calculate relatedness (r) based on a regression using the population sample of allele frequencies and assume HWE. Because no genetic structure was detected (see Results), the allele frequencies of the entire data set were used. Moreover, all estimators perform similarly when individuals are unrelated (Konovalov and Heg 2008). Pairwise relatedness estimators were calculated using the program KINGROUP version 2.0 (Konovalov et al. 2004). Mean r was calculated for all groups from different sampling locations and in different age cohorts. On average, r = 0 in unrelated individuals or when the relatedness of a group is random, r = 0.25 in half siblings, and r = 0.5 in full siblings. Pairwise relatedness was estimated for all recruits.
To determine the power of the 8 selected loci to detect full siblings, 17 000 null hypothesis dyads were generated with a Monte Carlo simulation, randomly reshuffling sample alleles at each locus following the method of Guo and Thompson (1992), also implemented in KINGROUP, which compares the observed and expected tables of genotypes to return a P value. A null hypothesis of no relationship was rejected if the P values were smaller than a prespecified α level that was estimated using the same pairwise false discovery rate as described above. Tests were repeated 10 times to gauge consistency. We also used KINGROUP to perform pairwise relatedness analysis of 12 pairs of simulated full siblings that had the same number of loci and allele frequencies as our sample data to calculate a type II error rate.

Molecular Diversity
Approximately 250 base pairs (bps) of the mtCR were resolved and these revealed elevated genetic diversity at this locus (Table 1). Overall, there were 164 variable sites, of which 17 were singletons, 53 were invariable sites, and 24 small, 1-to 2-bp indels. MtCR genetic diversity of recently recruited N. unicornis (h = 0.99, π = 0.085) was similar to that of the adult population on Guam (h = 0.99, π = 0.083). Genetic diversity at all recruitment sites was comparable and was in the range of h = 0.99-1.00 and π = 0.079-0.092. In total, 373 haplotypes were observed and 74.3% of all individuals had unique haplotypes. There was no significant genetic structure observed in the mtCR among sites, cohorts, or age classes after correction for a false discovery rate (Table 2). Overall, the data from the mtCR suggested that recruitment events of N. unicornis were unstructured in space and time.
Of the 12 microsatellite loci, only 2 had more than 5% and none more than 7% of missing data. Four of the 12 microsatellites (Nuni 02, 04, 09, and 10) had significant departures from HWE and were excluded from subsequent analyses (see Supplementary Material online). Out of the 28 pairwise tests of LD, only 2 were significant (loci Nuni 01 and Nuni 11, P = 0.03; Nuni 07 and Nuni 12, P = 0.04). However, these were based only on a single adult population. LD was also tested for recruits at each of the 14 sites, where 6 out of 392 tests were significant but never for the same locus pair. A relative lack of LD in recruits from the same site is inconsistent with kinship aggregations. According to MICROCHECKER, loci Nuni 02 and Nuni 09 were affected by null alleles, which may explain their departure from HWE expectations. At any rate, null alleles are expected to exaggerate genetic differentiation (Chapuis and Estoup 2007), which is clearly not the case here.
As with the mtCR, genetic diversity in the nuclear microsatellites was similar in both recruit and adult populations ( Table 3). The average inbreeding coefficient was low for most groupings (F IS = 0.012-0.066) and for cohorts at some collection sites, it was negative, indicating a lack of inbreeding and, by extension, a lack of relatedness within age cohorts from the same site. In total, 40 private alleles were observed. Private alleles may indicate that the genetic diversity is even greater than what our sample would suggest. Alternatively, rare alleles among recruits might be the result of genetic drift or possibly long-distance dispersal (Slatkin 1985). Pairwise Jost D values across 8 microsatellite loci revealed an absence of structuring among sites, cohorts, and age classes after correction for false discovery rate (Table 2). All hierarchal analyses of population differentiation produced nonsignificant values regardless of how populations were arranged (not shown).

Estimators of Relatedness
Pairwise relatedness between individuals and the mean relatedness coefficient for each group of cohorts is given in Figures 2 and 3. Overall, recruits had a mean relatedness of r = 0.0, which is expected under conditions of random mating. The mean for adults was r = −0.006. Mean recruit relatedness within each site ranged between r = −0.018 and −0.125. Notwithstanding a negative mean relatedness, some pairs from the same sites had high pairwise relatedness coefficients. For example, 1 pair, collected at Cocos, had a relatedness coefficient of r = 0.64 but did not possess the same

Table 2
Pairwise Fst values from mtDNA control region sequences (below the diagonal) and pairwise Jost D from 8 microsatellite loci (above the diagonal) for adult and recruit Naso unicornis individuals from Guam and Saipan, numbered as in Table 1 1 mtCR haplotype. For the purposes of this study, putative sibling pairs were defined only as those individuals that shared a mtCR haplotype (recruit mtCR haplotypes: N h = 295) and had a pairwise relatedness of r > 0.18, in at least one of the relatedness estimators. Although it is possible for true siblings to have a relatedness coefficient less than r = 0.18, we have arbitrarily chosen this cutoff value to be exclusive. Based on the above criteria, 12 putative sibling pairs were identified among sampled recruits (Table 4). Only 2 pairs, both from Asan, came from the same site but members of one of these pairs were sampled more than a year apart and so were not from the same cohort. All other putative sibling pairs came from different sites. Seven out of 12 putative sibling relationships were between-year classes. Two of the putative sibling pairs were composed of a Guam recruit and a Saipan recruit, separated by a distance of about 250 km. One of the Guam-Saipan pairs shared a relatedness of r = 0.29-0.62, enough to suspect a full-sibling relationship. This pair   was also the only putative sibling pair for which there was enough power to reject a null hypothesis of no relationship. This pair also shared an allele in 3 out of the 4 discarded loci that were not in HWE. In total, 271 pairs, about 0.3% of all comparisons, had P values small enough to reject a null hypothesis of no relationship (P < 0.0013). Yet, because all but one of these pairs did not possess the same mtCR haplotype, most of these were treated as unrelated individuals in this study. A small percentage of false positives are expected in any dataset (Weir et al. 2006), though some may be paternal half siblings (see Discussion). The type II error rate for our study, based on simulated full siblings, is estimated to be more than 20%. Therefore, although this study only claims enough power to identify 1 sibling pair, others in our sample data may have gone undetected.

DAPC Plots
Scatter plots of genetic variation in multivariate space revealed considerable overlap between Guam adults, recruits from 2008 and 2009 on Guam, and recruits from Saipan ( Figure 4). Each group occupies all quadrants of the plot and possesses individuals that lie well beyond the 95% inertia ellipse of the group. DAPC also failed to discriminate between groups of recruits when they were segregated according to sample site ( Figure 5), with the exception of Pago 2008 and Pago 2009, which do not overlap and mostly occupy separate quadrants. Otherwise, genetic variation in the microsatellite loci appears to be homogeneous in space and time.
The k-means algorithm identified 2 genetic clusters. Both clusters were composed of individuals from all age and year classes and included recruits from Saipan. When subjected to DAPC, these 2 genetic clusters were clearly differentiated by a single discriminant function (see Supplementary Material online). DAPC can sometimes create artifactual clusters (Jombart et al. 2010); therefore, the clusters may not be a biologically significant pattern. However, if the 2 clusters are a true pattern in the data, they are a signal that has arisen beyond the spatial and temporal scale of this study.

Discussion
Cohorts of settling marine organisms often exhibit genetic heterogeneity temporally, between settlement pulses, and spatially at microgeographic scales (Moberg and Burton 2000;Hellberg et al. 2002;Selkoe et al. 2006;Hepburn et al. 2009;Lee 2009;and others). From a metapopulation perspective, stochastic genetic variation, temporal heterogeneity, microgeographic structure, and sweepstakes reproduction suggest that migrant exchange between isolated adult populations is likely to be irregular, unpredictable, and potentially limiting, in terms of its demographic impact (Pringle and Wares 2007;Siegel et al. 2008;Cowen and Sponaugle 2009). Furthermore, increasing the time spent during the pelagic larval phase is expected to increase the irregularity of population connectivity on ecological time frames (Siegel et al. 2008). Hence, in a coral reef fish with a lengthy pelagic larval duration, such as in the case of N. unicornis (approximately 71 days), dramatic genetic variance between recruits sampled in time and space might be expected. To the contrary, however, we show that recruiting N. unicornis individuals display minimal temporal and spatial genetic structuring within the Marianas archipelago. It seems unlikely that passivity on the part of N. unicornis larvae could be responsible for the observed mixing, as these are some of the largest and most well developed of all reef fish larvae and have all the biomechanical specializations associated with high-performance swimming (Fisher and Hogan 2007). Therefore, some aspect of the larval biology of N. unicornis, such as behavior, must resist the natural forces that lead to temporal and spatial heterogeneity in other marine organisms.
Unlike the former study of , recruits were not sampled in crest nets on the night larvae entered the juvenile habitat and were instead collected days to weeks after settlement, which means that postsettlement mortality and genetic drift would have potentially changed alleles from their larval pool frequencies. Therefore, the genetic patterns presented here are not of larval pools per se. However, sampling genetic variation from postsettlement juveniles is arguably a better representation of the genetic contribution of each cohort to the adult population. Furthermore, after initial sampling, there was no observed change in genetic structure at settlement sites through time, as has been reported in other marine fishes Vigliola et al. 2007), and the earliest recruit samples of small recently settled individuals were just as genetically diverse as later collections containing larger individuals, as seen in Figure 3. Therefore, although understanding the genetic composition of larval pools has biological importance, how these pools converge upon settlement is more relevant to metapopulation considerations.
Considering that mixing and mortality after settlement are important demographic processes (Hixon and Webster 2002;Doherty et al. 2004), sweepstakes effects that may occur prior to settlement appear to be irrelevant to the genetic composition of the adult N. unicornis populations in the Marianas (see Supplementary Material online). Selkoe et al. (2006) came to a similar conclusion, stating that the temporal scale of heterogeneity in kelp bass (Paralabrax clathratus) was so fine that it had "no lasting effects." Other studies of marine fish have also shown sweepstakes effects to be minimal (Gilbert-Horvath et al. 2006;Priest et al. 2012). There is now increasing evidence that in many marine fishes, evidence of reproductive sweepstakes does not preclude high levels of metapopulation mixing (see Christie et al. 2010a).

Pairwise Kinship Relatedness
Microsatellite loci are valued for their hypervariability and fast mutation rates, which make them useful for investigating ecologically pertinent genetic patterns (Manel et al. 2005;Waples and Gaggiotti 2006). Nevertheless, the most polymorphic marker used in this study was the mtCR, which had 10.9× more variants than the most polymorphic microsatellite locus. Elevated polymorphism in the mtCR presupposes saturation at nucleotide substitution sites, possibly making it unsuitable for detecting patterns in evolutionary time. Yet, for ecological inferences, such as pairwise relatedness, polymorphism of the mtCR presents an opportune companion to microsatellite loci. Hence, kinship in this study was not inferred unless both individuals shared the same mtCR haplotype, which necessarily excludes all paternal half siblings but is otherwise an appropriate precaution.
Out of 74 097 pairwise comparisons between recruits, only 12 met the predetermined criteria for putative sibship ( Table 4). The chance that most of these pairs were true full siblings was low, particularly because only 8 loci do not have enough power to reject the null hypothesis of no relationship in most cases. One pair, however, received a relatedness score high enough to suspect full sibship and was statistically significant under Queller and Goodnight's relatedness estimator. Further, the mitochondrial haplotype shared by this pair had a relative frequency of only 0.0067 among recruits. Simplistically, because the relative frequency is also the probability that any 2 recruits will share that haplotype by chance, multiplying the haplotype frequency by the relatedness P value gives a rough approximation of the probability of a type I error (P = 0.0000038, about 1 in 264 000). Additionally, these 2 putative siblings shared alleles at 3 out of the 4 microsatellite loci that did not conform to HWE.
Overall, it appears plausible that these 2 individuals are at least maternal half siblings, although, it is unlikely that they were spawned at the same time because they were collected more than 3 months apart.
This single putative sibling pair that showed statistically significant results was collected from 2 sites, on 2 different islands, 250 km apart, further than any other kinship pair (including parent-offspring pairs) ever recorded in a coral reef fish. Other relatedness studies have shown that sibling reef fishes can settle at different sites Saenz-Agudelo et al. 2011), but until now, the scale of such observations has only been tens of kilometers. Here, we provide direct evidence that a single female can export offspring to multiple locations across an archipelago the size of the Marianas in the same year. Considering the limited sample sizes and limited number of loci used in our study, it might seem remarkable that we were able to detect any sibling pairs at all. However, our simulations indicate that the type II error rate among siblings may be greater than 20%. Thus, conceivably, a larger study with more loci would detect more siblings than we have, perhaps on a greater spatial scale also.
Correctly interpreting the results of sibling relatedness requires some knowledge of the reproductive biology of N. unicornis, which is not well known. Nevertheless, 2 insights are worth noting. 1) Mature males may form a mating hierarchy, as they sometimes aggregate in large numbers and afterward display wounds, apparently inflicted by the caudal spines of rivals (J.B.H., personal observations). If so, dominant males may have more mating opportunities than subordinate ones, leading to reproductive skew and possibly an abundance of paternal half siblings, which might explain the presence of recruits with high levels of relatedness that did not share a mtCR haplotype. 2) Although N. unicornis is known to aggregate, there is no evidence that mass spawning occurs. The only documented account of N. unicornis spawning is by Lieske and Myers (1994), who report that pair spawning was observed on the periphery of aggregations. Mated pairs of N. unicornis are also occasionally observed outside of aggregations (J.B.H., personal observations).
Some coral reef fishes mass spawn in large aggregations, where gametes from males and females are released with little regard for mate selection (Domeier and Colin 1997). This mingling of gametes is undoubtedly a source of population mixing for some species, but because N. unicornis does not appear to spawn this way, mixing must occur after the initial spawning.  reported that larval kinship groups were only evident after larvae were sorted by otolith age; the mean genetic relatedness of recruits settling on the same night was otherwise not different from random expectations. Therefore, substantial mixing must take place either during the pelagic larval phase or in the short time period when larvae relocate from the open ocean to the settlement habitat. The results of the present study may indicate that mixing primarily occurs before relocation to the settlement habitat, as the mean pairwise relatedness of recruits as a whole was much greater than the mean pairwise relatedness observed at any recruitment site (Figures 2 and 3). To the extent that the sampled recruits at different sites comprise nonindependent observations, this pattern suggests that the presettlement larval pool as a whole is a random assemblage of genotypes and that local postlarval recruit populations are, at times, less related on average than would be expected at random (Figure 3).
It is important to note that the large recruitment pulses sampled in this study are a feature of the population biology of N. unicornis not found in most reef fishes. Further, N. unicornis does not always settle in this fashion. To what extent our sampling of large recruitment pulses has affected our data is not known. However, our results suggest that pulses represent an admixture of larval propagules resembling a migrant pool (sensu Slatkin 1977).

Linking Demographic Processes with Broad-Scale Patterns of Connectivity
There is a perceived discrepancy between small-scale demographic patterns and broad-scale connectivity patterns in coral reef fishes and other pelagic dispersing marine organisms. It has been expressed that genetic connectivity studies overestimate, whereas demographic studies underestimate, migrant exchange (Hellberg et al. 2002). The underlying causes of this discrepancy have been labeled as 1) evolutionary versus ecological time-scale discordance, caused by the substitution rates of molecular markers (Hellberg 2007), 2) confounding factors such as large effective population size (Palumbi 2003, Hedgecock et al. 2007, and 3) biological differences in the target species of each type of investigation (e.g. reef fish demography studies are biased toward damselfishes). Here, we present a demographic study that does not downplay migrant exchange and is supportive of broad-scale phylogeographic patterns. Horne et al. (2008) reported that N. unicornis lacks genetic population structure across the tropical Indo-Pacific region, as far to the east as French Polynesia, and found shared haplotypes of the hypervariable mtCR at opposite ends of the geographic species range. The present study found a high degree of mixing among larval propagules, expected in a high-dispersal species, and even identified a sibling pair separated by 250 km. Therefore, investigations of N. unicornis at all spatial scales are in harmony.
All things considered, the upper limit of migrant exchange in N. unicornis appears to be beyond the scale of the Marianas islands alone. Computer simulations of population connectivity in some highly dispersive marine organisms suggest a bimodal dispersal kernel, where approximately 60% of recruits settle at distances less than 450 km from where they were spawned and approximately 22% settle over 1000 km away, with relatively few at intermediate distances (Butler et al. 2011). A bimodal dispersal kernel for the N. unicornis population on Guam seems not only plausible but also likely because shallow-water habitat is scarce at distances between 300 and 1000 km. Hypothetically, dispersal in N. unicornis across distances just a little more than 1000 km puts many areas in the range of Guam: Iwo Jima and the Ogasawara islands to the north, Yap and Palau to the southwest, and the Caroline Islands and Pohnpei to the southeast. Recruits sampled in this survey could have come from a large number of areas (see Mora et al. 2012), especially those lying to the east from which the North Equatorial Current flows (see Eble et al. 2011).

Conclusions
The bluespine unicornfish is both a useful species for investigating recruitment in coral reef fishes and an exploited fishery species (Rhodes et al. 2008), which plays a functionally important role on coral reefs as a consumer of coralsuppressing macroalgae (Hoey and Bellwood 2009). Spatially explicit conservation strategies, such as marine protected areas, are thought to be most effective when large numbers of demographic connections exist between subpopulations (Hastings and Botsford 2006). Such demographic connections are likely to exist for N. unicornis in the Marianas and other similar-sized archipelagos, considering that parents may export offspring to multiple destinations separated by approximately 250 km. Therefore, for N. unicornis and other highly dispersive reef fish species, networks of marine protected areas across hundreds of kilometers or more may be effective conservation and fisheries management tools.
Funding US Fish and Wildlife's Federal Assistance in Sportfish Restoration Program to the Guam Department of Agriculture; James Cook University.