Abstract

Integration of ecological and genetic data to study patterns of biological connectivity can aid in ecosystem-based management. Here we investigated connectivity of the Hawaiian grouper Epinephelus quernus, a species of management concern within the Main Hawaiian Islands (MHI), by comparing genetic analyses with simulated larval dispersal patterns across the species range in the Hawaiian Archipelago and Johnston Atoll. Larval simulations revealed higher dispersal from the MHI to the Northwestern Hawaiian Islands (NWHI) than in the opposite direction and evidence for a dispersal corridor between Johnston and the middle of the Hawaiian Archipelago. Genetic analyses using mitochondrial DNA (mtDNA) control region sequences and microsatellites revealed relatively high connectivity across the Hawaiian Archipelago, with the exception of genetically distinct populations and higher mtDNA diversity in the mid-Archipelago. These analyses support the preservation of the mid-archipelago as a source of genetic diversity and a region of connectivity with locations outside the Hawaiian Archipelago. Additionally, our evidence for directional dispersal away from the MHI lends caution to any management decisions that would rely on the NWHI replenishing depleted MHI stocks.

1. Introduction

Growing numbers of studies are integrating ecological and genetic data to investigate population connectivity [13]. These techniques have the potential to facilitate ecosystem-based management by increasing our understanding of the interactions between organisms and their environment. For marine species, these techniques may be particularly valuable because population genetic studies of these species are often thought to suffer from low power and therefore to be of limited value for management [46]. Studies of marine species commonly reveal weak population genetic structure that is described as “chaotic” due to its apparent unpredictability over space and time [710]. However, this weak structure may result from factors other than high genetic connectivity. For example, the presence of large population sizes in the marine environment may prevent genetic divergence due to the slow process by which large populations achieve migration-drift equilibrium [4, 11]. Additionally, violations of the assumptions of commonly used population genetic statistical models are thought to be particularly strong for marine populations; many of these populations have spatial and temporal heterogeneity in abundance and settlement rates which violate the assumptions of constant population sizes and migration rates typical of many population genetic models [5, 12, 13]. These violations are thought to result from the stochastic nature of oceanographic currents as well as cryptic spatial habitat heterogeneity in the sea [14].

Recent studies integrating environmental data with genetic data have proven useful in revealing the ecological drivers underlying the weak and apparently chaotic genetic structure of marine species. For example, genetic analyses indicated weak structure for the subtidal whelk Kelletia kelletii across the Santa Barbara Channel of California, but integration of genetic analyses with larval dispersal simulations based on ocean currents revealed that nearly 50 percent of this genetic structure could be accounted for by ocean currents [15]. Additionally, a comparison of genetic and environmental data for K. kelletii and two other marine species—kelp bass Paralabrax clathratus and California spiny lobster Panulirus interruptus—in the Southern California Bight indicated that kelp bed size was an important predictor of genetic structure for all species, with ocean circulation patterns also important for the kelp bass [14]. For each of these studies, the integration of environmental and genetic data provided a greater understanding of the ecological basis of population structure than could be gained from either data source alone. Ultimately these analyses provided a stronger basis from which effective management decisions could be made.

The Hawaiian grouper (Epinephelus quernus, Serranidae; recently revised to Hyporthodus quernus, Epinephelidae [16, 17]) is a species of management concern within the Hawaiian Archipelago. Endemic to the Hawaiian Archipelago and Johnston Atoll (Figure 1 [18]), E. quernus is part of a commercially important bottomfish species complex considered to have experienced long-term decline in the Main Hawaiian Islands (MHI) since the 1950s–1960s due to fishing pressures [19]. In the Northwestern Hawaiian Islands (NWHI), a phasing out the bottomfish fishery began in 2006 and was completed in 2010 as a result of the designation of this region as the Papahānaumokuākea Marine National Monument. Establishment of the Monument led to speculation as to whether dispersal from the protected NWHI to the over-fished MHI could replenish depleted MHI stocks. However, several life history traits of E. quernus may reduce the likelihood of high connectivity between the NWHI and MHI. Adults appear to be territorial, are found in relatively small groups, swim close to the bottom, have strong habitat preferences for rocky substrate within a depth range of 5–380 m (C.D.K., pers. obs.), and are protogynous [20, 21]. These observations suggest adults may form small stable spawning groups similar to other territorial protogynists. Additionally, the limited biogeographic range of E. quernus suggests that this species has overall low dispersal capabilities, although the correlation between dispersal ability and range size is not strong [22]. However in support of the NWHI as a possible source for MHI stocks, the relatively long pelagic larval dispersal (PLD) period of this species may promote connectivity facilitated by larval dispersal between islands. The PLD for this species is estimated at 35–45 days based on an analysis of otoliths from juveniles collected in the NWHI (R. Nichols and E. DeMartini pers. comm.); although this PLD is within the lower range of medium-to large-bodied groupers (range 30–80 days: [23, 24]), it may be long enough to promote connectivity across large geographic distances. The only data available to address this question comes from population genetic analyses using mitochondrial DNA (mtDNA) control region sequences for E. quernus samples collected throughout the Hawaiian Archipelago; these analyses found low or nonsignificant genetic divergence between most locations, with the exception of a few genetically distinct locations near the middle of the archipelago [25].

To investigate connectivity across the range of E. quernus using both ecological and genetic data, we compared simulated larval dispersal patterns based on oceanographic data with population genetic data. Our simulations utilized existing data on spawning behavior and PLD for E. quernus, as well as high-resolution data on ocean currents across the Hawaiian Archipelago and Johnston Atoll. Using data on interannual variability in current patterns over 17 years, we were able to simulate larval dispersal patterns across multiple generations, thus accounting for the cumulative effect of larval migration between intermediary sites over multiple generations [15, 26]. We then compared simulated larval dispersal patterns with the results of genetic structure analyses conducted for E. quernus samples collected across the Hawaiian Archipelago. For these analyses, we utilized previously published results from a study of mtDNA control region genetic structure [25], as well as new analyses using microsatellite loci for the same dataset as used in the mtDNA study. We then compared the results of larval dispersal simulations and population genetic analyses to make ecosystem-based management recommendations for E. quernus within the Hawaiian Archipelago.

2. Methods

2.1. Larval  Simulations

Larval transport simulations were performed using a version of a biased random-walk model as described by Polovina et al. [27]. The following modifications to this approach were made: firstly, surface currents used in this study were five-day interval (from January 1993 through December 2009) latitude/longitude resolution gridded flow fields from the Ocean Surface Current Analyses—Real Time (OSCAR) project at NOAA Earth and Space Research, henceforth referred to as OSCAR currents. The bounds of currents used in this study were from 170°E–140°W, 0°–50°N. These surface flow fields are a composite of altimetry-derived geostrophic components and satellite-derived wind components, tuned to 15-m depth drogue trajectories, thus capturing both the large-scale geostrophic motion as well as the surface, wind-driven Ekman transport. This approach is documented in Lagerloef et al. [28]. The OSCAR circulation data was used because it provided a large number of years ( ) to incorporate interannual and decadal oceanographic variability and allowed the use of as many years of circulation data as possible in a multigenerational simulation model. Secondly, simulations were conducted for both the lower and upper bounds of the estimated PLD for E. quernus (35–45 days, R. Nichols and E. DeMartini pers. comm.). These minimum and maximum PLD values were chosen for analysis due to evidence that these measures are better correlates of than is the average PLD for many marine species [29]. Thirdly, for spawning output, 100,000 larvae were released per year for the 17 years of available OSCAR data, ~85% of which were released during the spawning season of E. quernus from February through June, with a peak centered in March. A daily spawning output function was fit with nonlinear regression using a gamma probability distribution coupled with constant (representing both the spawning peak and low level of background spawning activity throughout the year); this was done to emulate the observed spawning patterns of E. quernus (R. Nichols and E. DeMartini pers. comm.). Fourthly, release and settlement habitat along the Hawaiian Archipelago was defined using a high-resolution bathymetric database [30]. A total of 959 habitat pixels were grouped to the nearest island, atoll, or bank, resulting in 26 locations in the Hawaiian Archipelago.

The value for eddy diffusivity used in this study was 500 m2 sec−1, the same as that used by Polovina et al. [27]. This value was originally derived by qualitatively tuning the larval transport simulation results to larval lobster catch patterns from a series of midwater trawls. Part of this study involved re-evaluation of this value of eddy diffusivity since the current field used in this study (OSCAR, geostrophic and Ekman transport) is different than the purely altimetric geostrophic currents used by Polovina et al. [27]. Eddy diffusivity was estimated from the trajectories of subsurface drifter buoys (drogued at 15 m depth) which occurred in the spatial domain of this analysis ( drifter buoys, daily positions = 38921, year = 2006). Drifter buoys are routinely deployed by research vessels and vessels of opportunity as part of a large program to understand global circulation patterns [31]. The drifter buoy dataset used in this analysis is administered by the Global Drifter Program at NOAA, and the data are available at: http://www.aoml.noaa.gov/phod/dac/index.php. Daily drifter buoy positions were merged with the OSCAR data to match the nearest value of currents to the buoy location. The predicted position for the next day was calculated using these currents, and this was compared to the actual position. By assuming that the discrepancy between these two values represents eddy diffusivity, the tabulation of differences can be compared to a tabulation of predicted diffusive steps using the Lagrangian formulation in the transport model. Three values of eddy diffusivity were examined: 100 m2 sec−1, 500 m2 sec−1, and 1000 m2 sec−1.

Multigenerational simulations were conducted by bootstrapping from the 17 available years of forward dispersal probabilities generated from the biased random-walk model and iteratively dispersing propagules for 1000 generations. Simplistic population dynamics (fixed population size with uniform mortality schedule) were assumed with a leveling imposed after each generation, while maintaining and tracking the new source breakdown. This evolving source breakdown incorporated inputs from the new generation added to what was there from the previous generations. The distributions equilibrated after a few hundred generations, and an average result was calculated from the final 100 runs, that is, generation 901–1000. Each island was fully self-seeding at the start of the simulation, that is, at generation 0. This simulation was run five times each for “forward” and “backward” projections to generate settlement probabilities for each island. “Forward” models report the final locations of larvae spawned at a given island, as measured by the percent of successfully settled popagules released from the source site that settled at the receiving site in a multigenerational context. “Backward” models report the natal islands of larvae which recruited successfully to each island, as measured by the percent of successfully settled propagules at the receiving site that originated at the source site in a multigenerational context.

2.2. Genetic  Data  Collection  and  Quality  Control

A total of 302 individuals were sampled from ten sites from across the Hawaiian Archipelago (Figure 1). Many of the specimens, particularly those from Kaua‘i through Pearl and Hermes, were obtained with the help of commercial fishermen operating in Northwestern Hawaiian Islands. Researchers from the Hawai‘i Institute of Marine Biology who were conducting habitat surveys provided most of the specimens from the Main Hawaiian Islands, and additional fish were purchased at the Hilo fish auctions to supplement sample size from the island of Hawai‘i.

For most individuals, latitude-longitude coordinates were recorded; however, this specific information was not available for 29 out of 36 specimens from the island of Hawai‘i that were purchased from the fish auctions. In addition, precise latitude-longitude data were not made available by fishermen for three out of 44 samples from Nihoa, 12 out of 20 from Gardner, two out of 46 samples from Maro/N. Hampton, and nine out of 30 specimens from Pioneer/Lisianski. For these individuals, the only locality information provided was the banks off the coasts of particular islands.

DNA was extracted using Qiagen DNeasy tissue extraction kits, and each individual was amplified using PCR for ten microsatellite loci. Five loci were among those isolated from E. quernus [32]; and the remaining five were developed for other epinephelids in the genus Mycteroperca [33, Zatcoff & Chapman, unpublished data] (Table 1). Several individuals were chosen at random and sequenced for each of the ten loci to confirm presence and characteristics of the repeat motifs. Once confirmed, 10 μL PCR reactions were performed on all samples using annealing temperatures shown in Table 1. Microsatellite fragments were then visualized on a CEQ 2000XL capillary system (Beckman Coulter). Allele sizes were scored and recorded for each locus and each individual using the CEQ 8000 Genetic Analysis System software 8.0 (Beckman Coulter). Deviations from Hardy-Weinberg Equilibrium were tested for each sampling location individually and for the total sample for all loci using ARLEQUIN 3.11 [34]. Each locus was also tested for null alleles, large allele dropout, and scoring errors using MICROCHECKER 2.2.0.3 [35]. Evidence for selection was examined for each locus using LOSITAN [36].

2.3. Genetic Diversity and Population Structure

Microsatellite diversity was examined by calculating observed and expected heterozygosities for each microsatellite locus, and each sampling site using the program ARLEQUIN 3.11. Allele richness was calculated using FSTAT [37]. Population genetic structure was investigated by calculating pairwise values across sampling sites using ARLEQUIN 3.11. Significance of pairwise comparisons was determined using 16,000 random permutations.

Population genetic structure across the Hawaiian Archipelago was further investigated using mtDNA control region sequences gathered from a previous study which utilized the same sample set as the study presented here [25]. Genealogical relationships between mtDNA haplotypes across the Hawaiian Archipelago were investigated by creating median-joining networks of mtDNA haplotypes using NETWORK [38]. Results from genetic analyses were qualitatively compared with results from larval dispersal modeling. Genetic structure was predicted to occur in regions with high larval retention or regions with asymmetric directional larval dispersal.

3. Results

3.1. Eddy  Diffusivity

The tabulation of drifter buoy diffusive daily steps was plotted with predicted diffusive steps at eddy diffusivity values of 100 m2 sec−1, 500 m2 sec−1, and 1000 m2 sec−1. The results indicate that the value of 500 m2 sec−1  best characterizes the daily diffusive spatial steps of the drifter buoys (Figure 2). The value of 100 m2 sec−1  predicted a narrower distribution of diffusion than that observed from the drifter buoy trajectories, while the value of 1000 m2 sec−1  predicted a much wider distribution of diffusion.

3.2. Simulated  Larval  Dispersal

Modeling of larval connectivity revealed similar patterns for 35- and 45-day PLDs (Figures 3 and 4), although as expected the 35-day PLD showed overall lower levels of connectivity between locations. Within the Hawaiian Archipelago, the highest levels of connectivity generally occurred between adjacent locations. However, forward and backward models for both PLDs indicated very little dispersal from the NWHI to the MHI (lower left quadrant of each graph), in contrast to relatively high dispersal from the MHI to the NWHI (upper right quadrant of each graph). Additionally, most sites throughout the Hawaiian Archipelago exhibited particularly high dispersal toward Gardner, and to a lesser extent, toward Necker and Ni‘ihau (Figures 3 and 4). Retention of natal larvae was also high at Gardner (Figures 3(a) and 4(a)).

For the 45-day PLD model, both forward and backward models indicated some connectivity between Johnston and the Hawaiian Archipelago (rightmost column of Figure 4(a); bottom row of Figure 4(b)), although as expected, dispersal between these two locations was lower than local retention (Figure 4). Both forward and backward models indicated higher connectivity between Johnston and Gardner than between Johnston and any other location in the Hawaiian Archipelago (Figure 4). Connectivity was also fairly high between Johnston and Necker (Figure 4). For the 35-day PLD model, there was no evidence of connectivity between Johnston and the Hawaiian Archipelago (Figure 3).

3.3. Microsatellite  Quality  Control

Microsatellite locus CA-4 showed evidence of null alleles for four out of ten sampling sites ( ) and deviated from Hardy-Weinberg equilibrium for six out of ten sampling sites. Microsatellite locus Gag010 showed evidence of null alleles for three out of ten sampling sites and deviated from Hardy-Weinberg equilibrium for four out of ten sampling sites. Tests for selection also indicated that locus Gag010 had a 99 percent probability of being under positive selection. To investigate the influence of CA-4 and Gag010 on analysis results, all analyses were conducted both with and without these two loci. The results differed little between analyses, and we report here the results without CA-4 and Gag010. No other loci showed evidence of null alleles, selection, or deviation from Hardy-Weinberg equilibrium.

3.4. Genetic  Diversity  and  Population  Structure

Observed and expected heterozygosity values, average number of alleles, and allele richness for microsatellite loci were similar across the Hawaiian Archipelago (Table 2). For mtDNA control region sequences, genetic diversity was higher in the NWHI than in the MHI, with the highest value occurring at Gardner (Table 2).

Genetic structure analyses using microsatellite data revealed significant pairwise values for six out of 45 pairwise comparisons ( ; ) (Table 3). Each significant pairwise comparison included at least one sampling location in the mid-archipelago (Pioneer/Lisianski through Nihoa). For mtDNA analyses, 11 out of 45 pairwise comparisons were statistically significant ( ; ), with each of these significant comparisons also including at least one sampling location in the mid-archipelago (Table 3).

Eight shortest mtDNA haplotype networks were found, but each network was similar in structure. One of the shortest networks was randomly chosen to report here (Figure 5). The network revealed that many haplotypes were shared across the archipelago, and most haplotypes were only one or two base pairs different from other haplotypes. The most common haplotype was found at all sampling locations except Gardner, and this haplotype was only one base pair different from many other haplotypes. However, several haplotypes or haplotype groups were at least six to eight base pairs different from any other haplotypes; most of these divergent haplotypes were found in the NWHI, although some were also found at the island of Hawai‘i (also known as the Big Island) and Maui. Notably, four of the most divergent haplotypes occurred at Gardner despite the relatively small sample size at this location ( ).

4. Discussion

Both genetic analyses and larval modeling revealed relatively high connectivity across the Hawaiian Archipelago for E. quernus. Limited population structure appears to be a common feature for most marine fish species studied in the Hawaiian Archipelago thus far. For example, no genetic structure was found throughout the Hawaiian Archipelago for bigscale soldierfish Myripristis berndti [39] or surgeonfish Acanthurus nigrofuscus [40], and genetic divergence was found for only one sampling location (Kure Atoll) within the archipelago for damselfish Stegastes fasciolatus [41]. Other species have exhibited low but significant genetic structure across the archipelago, including two surgeonfish: Ctenochaetus strigosus (average pairwise ) and Zebrasoma flavescens (average pairwise ) [40]. As a counter example, however, the damselfish Dascyllus albisella had strong population genetic structure across the archipelago, with pairwise values ranging as high as 0.718 [41]. The strong structure for D. albisella was proposed to result from the interaction of species-specific biological factors (spawning seasonality and habitat preferences) with environmental factors (seasonal ocean current patterns) [41].

Given evidence for a benthic adult lifestyle and a 35–45 day pelagic larval dispersal period for E. quernus, connectivity across the archipelago for this species is likely driven by larval dispersal rather than adult migration. Therefore oceanographic current patterns are expected to have a strong influence on dispersal patterns for this species [42]. Our larval dispersal simulations based on oceanographic currents across the Hawaiian Archipelago and Johnston Atoll provide evidence that dispersal within the archipelago may be dominated by dispersal from the MHI toward the NWHI and that dispersal in the opposite direction is probably low. Evidence for directional dispersal from the MHI to the NWHI has been observed for other marine species as well. For example, a previous study simulating larval dispersal for corals with a PLD of 60 days using a two-dimensional Eulerian advection-diffusion model also found predominantly northwesterly connectivity across the Hawaiian Archipelago [42, 43]. Additionally, population genetic analyses indicated higher gene flow from the MHI to the NWHI than in the opposite direction for both limpets (Cellana spp. [44]) and sea cucumbers (Holothuria atra [45]). This directional dispersal provides evidence that the North Hawaiian Ridge Current (NHRC) exerts a strong influence on larval movements for E. quernus; the NHRC flows northwest along the MHI [46] and was also proposed to influence genetic structure in this species by Rivera et al. [25]. The evidence for directional dispersal revealed by our larval modeling analyses does not support the proposition that the protected NWHI could act as a source to replenish the depleted MHI bottomfish fishery.

Both genetic analyses and larval modeling also revealed unusual patterns of connectivity for two sites near the middle of the archipelago: Gardner and Necker. Population genetic analyses revealed genetic distinctions between these two sites versus many other sites in the archipelago, and the highest mtDNA diversity value occurred at Gardner. Examination of the mtDNA haplotype network revealed that the high diversity at Gardner resulted from the presence of a large proportion of highly divergent haplotypes and the absence of the haplotype found most often in the rest of the archipelago. The presence of divergent haplotypes at Gardner may result from recent connectivity with a source outside of the Hawaiian Archipelago; this source would most likely be Johnston Atoll, because this is the only location outside of the Hawaiian Archipelago where this species has been reported. Unfortunately, E. quernus samples were not available from Johnston Atoll for this study, and we could not evaluate this hypothesis. Alternatively, the divergent haplotypes at Gardner may be remnants from an ancestral population, and their presence may be evidence that the ancestral population for the Hawaiian Archipelago occurred at Gardner, as suggested by Rivera et al. [25]. Other studies [44, 45] have highlighted unusual patterns of genetic diversity in the central NWHI (Gardner, French Frigate Shoals and/or Necker) indicating that this finding is not a spurious result. Whether these divergent haplotypes are a result of recent connectivity or ancestral haplotypes, both of these scenarios provide evidence that the central NWHI, in the area of Gardner/French Frigate Shoals/Necker, may be an important source of genetic diversity for the Hawaiian Archipelago.

Larval modeling indicated that most sites in the archipelago had unusually high connectivity with Gardner and to a lesser extent Necker and Nihoa, with all models indicating dispersal from most sites toward these three locations. For the 45-day PLD model, both forward and backward models also revealed higher levels of connectivity between Johnston Atoll and Gardner/Necker than between Johnston Atoll and any other locations within the Hawaiian Archipelago; again in these cases, the level of connectivity was higher for Gardner than for Necker. These analyses provide evidence that the genetic divergence of Gardner and Necker may result from the accumulation of genetic material through connectivity with multiple sources across the Hawaiian Archipelago as well as with Johnston. A connection between Johnston and the middle of the Hawaiian Archipelago near Gardner has been proposed for other marine species as well based on biogeographic or genetic data, including corals within the genus Acropora [47], the crown-of-thorns seastar Acanthaster planci [48], and vermetid snails [49]. A previous study simulating larval dispersal patterns based on ocean current data also found evidence for a dispersal corridor between Johnston and the middle of the Hawaiian Archipelago near Gardner for species with a PLD greater than 40 days, facilitated by the subtropical countercurrent and the Hawaiian Lee countercurrent [50]. Similarly, a study simulating larval dispersal for corals with a PLD of 60 days also found higher connectivity between Johnston versus the mid-archipelago (Maro, Gardner, and Necker) than Johnston versus other sites within Hawai‘i [42, 43]. Our study provides further evidence that the middle of the Hawaiian Archipelago is a region of connectivity with locations outside of the archipelago.

5. Conclusions

Our comparative analyses of genetic and oceanographic data indicated fairly high connectivity across the Hawaiian Archipelago for E. quernus, driven primarily by directional dispersal from the MHI to the NWHI. The mid-archipelago around Gardner and Necker was found to harbor higher genetic diversity than the rest of the archipelago, possibly due to connectivity between this region and Johnston Atoll. Our evidence for low levels of directional dispersal from the NWHI to the MHI lends caution to any management efforts that would rely on the NWHI replenishing depleted stocks in the MHI. Additionally, our analyses highlight the importance of the central NWHI from Gardner to Necker and support preserving the mid-archipelago as a source of genetic diversity for this species. Overall these analyses illustrate the benefits that can be gained for ecosystem-based management by integrating both environmental and genetic data for studies of biological connectivity.

Acknowledgments

First authorship is shared between M. A. J. Rivera who collected all genetic data and conducted preliminary genetic analyses and K. R. Andrews who conducted additional genetic analyses in the context of oceanographic simulations. The authors thank V. Moriwake, A. Moriwake, B. Alexander, H. Richman, M. Chow, and E. G. Grau for their help in collecting specimens and technical support. They also thank John Reed, Greg Holzman, Silas Naig, Dane Johnson, and Guy Ohara who are commercial bottomfishermen of NWHI (now closed) and MHI. They thank Ryan Nichols, Edward DeMartini, and Robert Humphreys from the Life History Program at PIFSC for sharing unpublished information on E. quernus life history. The authors thank Crow White for tutoring them in the use of his Matlab code for estimating oceanographic distance. Funding was provided by the Hawai‘i Department of Land and Natural Resources (Division of Aquatic Resources), Western Pacific Regional Fishery Management Council, University of California, Berkeley, National Science Foundation (DEB#99-75287, OCE#04-54873, OCE#06-23678, OCE#09-29031), Office of National Marine Sanctuaries PMNM-HIMB partnership (MOA-2005-008/6882), and the University of Hawai‘i Sea Grant College Program. Genetic data was part of the PhD dissertation of Malia Rivera while at UC Berkeley. This is HIMB contribution no. 1411 and SOEST contribution no. 8026.